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PREFACE 


This report is an edited version of notes distributed at the Summer Work- 
shop on the UCLA General Circulation Model in June 1971. It presents the 
computational schemes of the UCLA model, along with the mathematical and 
physical principles on which these schemes are based. 

Included are the finite difference schemes for the governing fluid- 
dynamical equations, designed to maintain the important integral constraints 
and dispersion characteristics of the motion. 

Also given are the principles of parameterization of cumulus convection 
by an ensemble of identical clouds. (A newer parameterization of cumulus 
convection, by an ensemble of clouds with a spectral distribution of sizes, 
will be published in a subsequent report.) 

A model of the ground hydrology, involving the liquid, ice and snow 
states of water, is included. 

A short summary is given of the scheme for computing solar and infrared 
radiation transfers through clear and cloudy air. A more detailed description 
of the scheme, by A. Katayama, is in press as Technical Report No. 6. 

The research reported here was done with the support of the National 
Science Foundation, Atmospheric Sciences Section (Grants GA-1470 and 
GA 34306X) # and the National Aeronautics and Space Administration, Institute 
for Space Studies, Goddard Space Flight Center (Grant NGR 05-007-328). 
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I. GOVERNING EQUATIONS 
IN o-COORDINATE 

1 . a -coordinate. 


The vertical coordinate used in the model is a, defined by 


Ps “Pt 


( 1 . 1 ) 


where p T is the pressure at the upper boundary of the model atmosphere, taken 
as a constant, and p s is the pressure at the earth's surface, which is a function 
of the horizontal coordinates and time. It follows from (1. 1) that 


a = 0 at p-p, , (|.2) 

0 = 1 at p = p s , J 

so that the boundaries are coordinate surfaces. 

i 

I b 

We define tt by 

ir = p s - p T • (1.3) 

ir/g is the mass of the entire vertical column of the model atmosphere with unit 
horizontal cross section. From. (1.1) and (1.3), 



p = p T + UO . 


(1.4) 
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From (1.4), 


(0 


= s!e = 

dt 


ircr + <; 


dir 
dt ' 


where cr = da/dt. Since it is a function of the horizontal coordinates and time 
only. 


to - ira + (r( +\V *v)ir 


(1.5) 


where W is the horizontal velocity and v is the horizontal gradient operator. 

We must remember that ircr is not identical with to. 

The earth's surface is a material surface as well as a coordinate surface. 
The kinematical boundary condition there is simply a = 0. At v = 0 (p = p T ), we 
have ircr = u)p = p . We may assume that to^ -0 is an acceptable approximation, 
if p T is chosen properly. Then 

a = 0 at a = 0,1 . (1.6) 


2. The hydrostatic equation . 

When pressure is the vertical coordinate, we have 



where <$ = gz and a is the specific volume. From (1.4), 


(1.7) 


- 2 - 



1-3 


6p = ir 6cr , (1.8) 

where 6 denotes the differential under constant horizontal coordinates and 
time (and, therefore/ under constant ir). (1.7) and (1.8) give 

5$ = -tt a 6a . (1*9) 

The following alternative forms of the hydrostatic equation can be 
derived from (1.9) and will be useful: 


6($or) = (4>-tMra)6a , 

(1.10) 

6$ = - c p 0 6(p X ) / 

(1.11) 

6(c p T+4») = p x c P 69 , 

(1.12) 


where x = R/c P and 0 = T/p X . 

3. The equation of continuity. 

In the pressure coordinate system, the equation of continuity takes the 

form 

v p -\V+|^=°. (U3) 

We have the relation 

v p = v a +(v P a)-^- . (1.14) 

This relation is illustrated in the accompanying figure. 
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cr 

c r + a<t 


Let A be an arbitrary quantity. Then 


A a + A Q " A h Act 
A s As Act As 


The limit as As -* 0 gives 


a A . 9A 
v p A = v A + x~ Vpcr 

p CT OCT p 


From (1.1) and (1.3), we obtain 


V P cr = Vp(£^-) = 


- — Vir . 
IT 


Then (1.14) gives 


v p = v — Vir x— . 


ct ir 


3ct 


(M5) 


Using (1.15) for v P *\V and (1.5) and (1.8) for 3cju/3p, results in 

cr 3\Vi , 3 r • . / 3 




and finally 

+ v^ «(ir W ) + ^ (iro ) = 0 . (1.16) 

This form of the continuity equation, with the cr -coordinate, is similar to the 
continuity equation with the z -coord inate , 
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|^ + V • (p W) + ^ (pu)) - 0 , 

where p is density. Let AS be a horizontal area element. Then p Az AS is 
the mass of the volume element Az AS. Similarly/ (rr/gJAa AS is the mass of 
the "volume" element Act AS, in a-space. However, tt is a function of the 
horizontal coordinates and time only, whereas p may also be a function of the 
vertical coordinate. 

The equation of continuity (1.16) is used for two purposes: to find 

3ir/3t (which is 3p s /3t) and to find a. 

Integrating (1.16) with respect to a, from 0 to 1 , and using the boundary 
condition (1.6), we obtain 


3ir 

3t 


= - J V (tt W ) d a 

i 

= - V • J ir W d or . 


(U7) 

(1.17)' 


This is equivalent to the so-called "surface pressure tendency equation". 

Substituting 3ir/3t, obtained from (1.17), into (1.16), we get 3(ir<7)/3cr. 
Integration of this derivative with respect to a , from 0 (or 1) to a given cr, gives 
-rr cr at that a-level. 


4. The individual time derivative and its flux form. 


With the a -coordinate, the individual time derivative d/dt is 
expressed as 
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-jr = (;sjr) +\V*v + o~ . 

dt ot cr a acr 


(1.18) 


Let A be an arbitrary scalar. By using the continuity eq. (1.16), 
we obtain the flux forms 

ir F = (l r )r (irA) + V (irWA) + ^' (Tr ‘ TA) ' 


(1.19) 


and 


*A ¥ = (l")a (^ A2 ) + ^ a •(* A2 ) +£(w* A* ) • (1*19)' 


In (1 . 19)' , A may be a vector. 


5. The equation of motion. 

The pressure gradient force is given by - v P $. Applying (1.15) to $ , 
we obtain 

v P $ = v , (1.20) 

p <7 IT OCT 

and substituting from (1.9) into (1.20), 

V p $ = V a <l>+ aa v tt . (1.21) 

In the CT-coordinate system, the pressure gradient force consists of 
two terms, as shown by equation (1.21). Where the slope of the earth's surface 
is steep, the individual terms are large but are approximately in the opposite 
directions. In the particular case where v p <l>- 0, complete compensation occurs, 
as shown in the figure below. 
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2 

<r- surface 

l 




The horizontal component of the equation of motion is 



~+/lkx\V + V<i> + craVTr=/P , 
at o 


( 1 . 22 ) 


where {F is the horizontal frictional force and d w/dt is the horizontal acceleration. 
Note that 

ir(v $ + aavir) = V (ir$)-($-<Tccn)viT.. 
cr o 

Using (1.10), $ - a aw = 3($a)/3cr. This gives us another form of the equation 
of motion. 


*cF + x + v gj (**)?* = w 


F , ( 1 . 22 )' 


6. The equation of state . 


a = 


RT 


(1.23) 


7. The first law of thermodynamics . 


The specific entropy is c P 6m.Q + const. , where 0 = T/p X . The first 
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law of thermodynamics is 


or 


d „ _ Q 

j|. Cp Q/n 0 j / 


dQ 1 0 n 

dt c P T Q ' 


where Q is the heating rate per unit mass. 

The flux form which corresponds to (1 .24) 1 is 


5T<*8> + v (* we) +£ (we) - r- f Q 


The first law of thermodynamics can also be written as 


dt 


c p T = uja + Q , 


where 


uj = ira +a(g^- + W • v) ir , 


as was given by (1.5). 

The corresponding flux form is 


■gy (irCpT) + v.(ir\Vc P T) + a^( 1T< ^ c P T ) - IT 0)0 +TTQ 


3cr 


Using (1.5), (1 . 25) 1 may be rewritten as 


■|^(ttc p T) +V.(ir \Vc p T)+p X c P 0) - Trcrct(^-- + W *v)rr + ttQ 


da 


(1.24) 


(1.24)' 


(1.24)" 


(1.25) 


(1.25)' 


(1.26) 
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8. The water vapor equation . 

Let q be the water vapor mixing ratio. The continuity equation for 
water vapor is 

^ = -C + E , (1.27) 

where C is the rate of condensation and E is the rate of evaporation per unit 
mass of dry air. The corresponding flux form is 

■|j(irq) + Vj(ir\Vq) + ^ (iraq) =ir(-C+E) . (1.28) 


II. INTEGRAL PROPERTIES 


The following integral properties of the governing equations, or of 
selected terms in these equations, are useful in designing the vertical finite 
difference scheme. 


1 . Mass conservation . 
(1.17)' gives 



i 

V. J irWda . 

o 


(II. 1) 
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The area integral of (11.1) over the entire globe vanishes. This shows that the 
total mass of the model atmosphere is conserved. 

2. Vertically integrated horizontal pressure gradient force . 

With the z-coordinate, the horizontal pressure force per unit mass is 

V z p, and per unit volume it is -v 2 p. Integrating vertically, 

P 

00 p 03 -» 

- J* v z pdz = -|vj pdz + p s vz s ! , (11.2) 

z s z s 

where z s is the height of the earth's surface. 

A line integral of the tangential component of the first term on the right 
in (U. 2), taken along an arbitrary closed curve on the sphere, always vanishes; 
and only the second term can contribute to the line integral. Therefore (except 
for the possible effect of a surface stress), only when there is a non-horizontal 
boundary surface can there be a "spin-up" or "spin-down" of the atmosphere 
along the arbitrary curve. 



The accompnaying figure is a vertical cross section of the earth's 
topography in the plane of its slope. p s n is the pressure force normal to the 
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IT 

surface and the horizontal component of this force is p 3 cos (2 "<=>) = Ps sina. 

This is the horizontal force per unit area of the surface. Per unit horizontal area, 
we have p s sina/cosa = p s tana = p s | vz s [ . 

With the p-coordinate, the horizontal pressure gradient force per unit 

mass is -V p $. Its vertical integration, with respect to mass, is 

i Ps i Ps 

— f V p <S>dp = — F V P $dp -<i> s v p s (11*3) 

g J o g L J o 

1 Ps 

= — Tv r («i>-$ s )dp +p s v$ s ] . 

With the cr -coordinate, if we start with the form of the horizontal pressure 
gradient force given in (1 .22) 1 , we immediately obtain the equivalent relation 

~gJ* [ v a ( 1T $)“ Virjda = - ^ jvJir($-^)da+Trv^ s j . (11.4) 

In order to derive (11.4) from the form given by (1.22), we must use the 

relation 

$3 = r ($ - CTTT a) dcr , (11*5) 

d o 

which is obtained from (1.10). (II. 5) can be rewritten as 

J (<J>-$ s )dcr = J* airadcr . (II .5)* 

0 0 o 

When p T = 0, (II. 5)' is equivalent to the familiar relationship between the vertically 
integrated geopotential and internal energies. 
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3. The kinetic energy equation. 


From the equation of motion (1.22), we obtain 

ir \V 2 = -ir\V • $ + a a vtt ”1 + ir W • (P , 


(H.6) 


The left hand side of (11.6) can be written in the flux form (1.19)'. 

The rate of kinetic energy generation by the pressure gradient force is 


-it\V. V $ + a a 7 ir 
La 


-7 *(ir \V$) + $7. (it \v) - aira W *7ir 
a a 

-7 *(tt W $) - $ I 57 - + 5 — (ira)l-OTra \V*7ir 
a Lot da J 

9 0ir * 0 $ 

-7 .(irW^)- (iTd<3>) -$ x— +ira r oira\V*7ir 

a da dt da 

- 7 o ,.(Tr\V$)-^(TTa^-«?-aTra)|^ 

- ir Fc (|^ + W • 7i T ) + ira ja 

- 7^*(tt\V$)- ( ir o - $ + <!? a |^) - ir uoa . (11.7) 


Here (1.16), (1.9), (1.10) and (1.5) were used. The vertical integration of the 
last form of (1 1.7) is 

- — f Tr\V*r7<l>+aa7Trjda = - — 7 • [ ir W $ da - — <l?s Ir “ “ F ir uo a da . 

g J o La J g J 0 g 9t g J 0 


4. Conservation of total energy. 


From (11.6) , (11.7) and (1.25 ) 1 , we obtain 
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(Jr) w2 ) + v a '( 1T w|\v 2 )+ ^ (ttc^W 2 ) = -TrW'fv^'J' +aa v ir~j +ir W-F, (11.9) 

1 cr 

• (tt W <£>) + ["(ff + ira)$j = ir\vTv ff $ +cravirj -Trcua , (11.10) 

and 

^-(irCpT) +V ^ , (tt \Vc p T) +^- (ira c p T) = irQ +ira)a . (M.ll) 

Taking the sum of (11.9), (11.10) and (II. 1 1) , and integrating with 
respect to a, from 0 to 1 , gives 

+J ir(£w 2 + CpT)daJ + V« J^irV/d \V 2 +c p T +$)da 

= f n(\V'p+Q)dcr . (11*12) 

‘ J o 

The area integral of (11.12) over the entire globe makes the contribution of the 
divergence term vanish, and we have the conservation of total energy when 
F =0and Q =0. 

5. Integral constraints on 8 and 0 3 . 

Under an adiabatic process, we have d0/dt = 0. The corresponding 
flux form, given by (1.24)", is 

^(ttQ) +V a .(ir\V9)+ £(md) = 0 . (11.13) 

Integration of (11.13) with respect to a, from 0 to 1 , gives 

•Ir f ir0da + v • f ir \V 0 da = 0 . (11.14) 

Or * o o 

Because the second term in (11.14) vanishes when the area integral is taken over 
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the entire globe, we see that the global integral of 0, with respect to mass, 
is conserved. 

Because (d/dt)0 2 = 0 under an adiabatic process, we can similarly 

derive 

P ir0 2 da +v • [ ir\V0 2 dcr = 0 . (11.15) 

Again the second term vanishes when the area integral is taken over the entire 
globe so that the global integral of 0 2 , with respect to mass, is conserved. 

Under an adiabatic process, the frequency distribution of potential 
temperature does not change with time. The integral constraints are not suf- 
ficient to maintain the frequency distribution of 0, but they effectively maintain 
its variance as well as its mean. 

6. Introduction to vertical differencing . 

The integral properties discussed above will be used for the design of 
the vertical finite difference scheme that is presented in the next chapter. 

The solutions obtained with any convergent scheme will satisfy these 
integral properties in the limit as the vertical grid size approaches zero. But 
the solutions obtained with these various schemes approach the true solution 
through different paths in a function space. Our aim is to find that scheme 
whose solution approaches the true solution through that path along which the 
finite difference analogs of the integral properties are maintained regardless of 
the grid size. 
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The principle that will be used here is similar to what is usually done 
when approximations are introduced into the governing equations. For example, 
when the hydrostatic equation is used as an approximation, we also use an 
approximate form of the horizontal component of the equation of motion, so 
that certain integral properties, such as energy conservation, are maintained. 

In doing that, however, the definition of energy is changed from its original 
definition. This modified energy approaches the true energy as the accuracy 
of the hydrostatic approximation increases. 

Many schemes which do not maintain the integral constraints on 0 and 
Q 2 were designed and tested. But none of these gave better results with long 
term integrations than the scheme, described in the next chapter, which does 
maintain these integral constraints. 
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III. VERTICAL DIFFERENCING, PART A 


1. The vertical index . 

The index k is used to identify a level. 

At levels with odd k, \V, T and q are carried. 

At levels with even k, d is carried. 

The upper boundary is k = 0 and the lower boundary is K+ 1 . 

0 /////////// / / / / / /-Z 4=0 „=0 

1 \y, T , q 

2 — a a = a 3 

3 \V,T,q 

k-2 \v, T , q 

k-1 — 6 CT=(J i c - 1 

k W,T,q 

k+1 — ° CT=ff k+i 

k+2 W,T,q 

K-2 W,T,q 

K-1 — d o= cr K _ x 

K \V,T,q 

K+1 ////-/-///// 7 /////// * = 0 a = 1 

We define ^ cr k= a k+x ” ' 0^*1) 

Then i/Acr^ = 1 , where £' is summation over odd k. We also define 

k =1 

*k = !K+i +cr k-i) • t 111 ' 1 )' 

The current UCLA 3-level model (December, 1971) uses p T = 100 mb, 

K = 5, Aa 1 = 3/9, A<t 3 = 4/9 and Aa 5 = 2/9. 


- 16 - 



ill-2 


2. The equation of continuity . 

We use the form 

|^+v(ir\V k ) + = 0 * ( ,IL2 ) 

S' (lll.2)Aa k gives 

k = l 

which is an analog of (1.17). 
mx k+r is given by 

w k+1 = [If + v(Tr\V k )]Aa k . (III. 4) 


lf = - S' V.(ir\V k )Aa, 

OT k =i 


(III. 3) 


3. Flux forms . 

For a variable A, carried at odd levels, flux forms analogous to (1.19) 
can be written as 

-^(uAj +v-(Tr\V k A k ) + ^—■(Tra k+1 A k+1 -iro' k _ 1 A k _ 1 ) . (III. 5) 

Because A is a variable carried at odd levels, we must somehow determine A 
at even levels (such as A k+1 , A^). 

From (III. 5) and (III. 2) we get 

"fir* + W.-vA. + ^-^A.^-A,) +d k _ i (A k -A k _ 1 )) ! . (III. 6) 

The expression in the bracket gives the form for dA/dt which is consistent with 
the flux form (III. 5). 
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If we require that nA dA/dt also be written in a flux form corresponding 
to (1.19)' 7 then A at even levels must be chosen properly. (III. 6) multiplied by 
A k gives 

''[^■4A?+W,-v|A| + 2 J-(d ktl (A K1 . 1 A l -A|)+a k . 1 (At-A k A k . I ))]. (III. 7) 
Using the equation of continuity/ (III. 7) can be rewritten as 


fj-faiA*) + V*(Tr\V k |A 2 ) +^- [ua k+1 (A lc+1 A k -|Af ) 

-ir<r k _ 1 (A k _ 1 A k -|A2) 


(III. 8) 


In order that (III. 8) be in flux form, A k+ a A k -^ Af must be i?A 2 at level k+1 
and A k _ 1 A k -^ A^ must be ^A 2 at level k-1 . Therefore,A k _ 1 A k -^A| , with k 

replaced by k+2 r must be equal to A k+ 1 A k - §A^ . Then 

A A -iA 2 =A A -iAs 

ft k+l tt l!+3 2 M k+ 3 M k + 1 ^k 2^k ' 

or 

Ak+i (A k + 2 “ A k ) _ 1? (A k + 2 “ A k ) • 


Because A k+3 - A k is generally not zero, we must have 

Ak+1 ~ +i ^k + 2 ) ' 


4. The acceleration term . 

Following (III. 6), we write the acceleration term as 

^ "g^ + (^ic‘v)\V k f^k+i(^k + i“ W k ) + 0 ' k _ 1 (\v k - \V k _ x . ( 1 1 1 p 1 0) 

To have a flux form for \v k -(d\v/dt) k , we choose 

w Jc+1 =i(w lc +\v k+3 ) . (HMD 
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This guarantees the conservation of total kinetic energy, insofar as vertical 
advection is concerned. The finite difference expression for the kinetic 
energy in a vertical column is 


— S' §W,f Ao k . (111.12) 

g k=i 

K 

S' is the summation for odd k. 

k =1 


5. The pressure gradient force . 


We wish to maintain the property of the vertically integrated pressure 

gradient force discussed in Chapter II, Sec. 2. For this purpose, it is convenient 

3 

to start from the form given by (1.22)' . We write v a ( ir ^ > )" fo (^u)vir as 

| A A 

• (NI-13) 

The symbol * is a reminder that the variable is at an even level. The analog 
to (11.4) is 


-1 S' (111.13) At,=-£[V( 

g k=i y “ 1C=1 ' 5 J 


In this way the integral property is maintained. 


(111.13) is equivalent to 


(111.14) 


Let 


* ‘ ^ & +i a k+i “ vir ]. (111.15) 

»(oa) k s# k - XT" - ^k-i^k-i) • (111.16) 

LMJ ^ 
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This is an analog to (1.10). However, at this stage (111.16) is only the definition 
of the symbol (<Ta) k , rather than the hydrostatic equation, because the dependency 
of (aa) k on temperature, pressure and a is not yet specified. Using (III. 16), the 
quantity inside the bracket of (111.15) can be written as 

V$ k + (aa) k Vtr , (111.17) 

which is an analog of V +oraVir at level k. 

6. Kinetic energy generation . 

To obtain the kinetic energy generation in finite difference form, we 
follow the procedure used in deriving (11.7). 

- * W * '[v$ k +(cra) k virj = 

0 | 

= - V-(Tr\V k $ k )-$ k [^+^-(TO k+ 1 -TO k _ 1 )J--ir((ra) it \V k ' Vir 
= - V.(Ti\V k $ k )-^- (Tra k+ 1 $ k+ 1 -ira k _ 1 $ k _ 1 )- $ k ^ 

1 p * ~ i 

+ + “^lc -i> | -ir(aa) k \V k -Vir 

1 a a g^. 

= - V.(ir\V k <l» k )-^-(TO k+1 <|) k+1 -TO k _ 1 $ k _ x ) - ($ k -ir((m) k )g^- 

p3 lf A 

-ir[(aa) k (^- + \V k -v)Tr- {ira k+1 ($ k+1 -$ k ) + Tia k _ 1 ($ k - 5> k _ x )} j 

= - V-(Tt\V k * k ) - [(Tta k+1 +a k +1 |f) $ k+1 - (w k _ 1 +a k _ 1 |f ) j 

- Tr(uua) k . (111.18) 

Here (wa) k is defined by 
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(iua) k = (cr a) k (^" + W k *v)tt- ^^{^k+i^k +i“$k) +ir ^ic-i(^ , k“^k- 1 )} • (HI. 19) 

At this stage, (111.19) is the definition of the symbol (cua) k . 

From a finite difference scheme for the first law of thermodynamics, we 
will determine a form for ( a>a) K . Then, by comparing it with (111.19), we will 
determine a form of (<ro) k and then a finite difference expression for the hydrostatic 
equation. 


7. The first law of thermodynamics. 


As an analog to (1.24)", we use the form given by (III. 5) with A - 0 
+ V.(Tr\V k 0 k ) + ^-(TWk+i0k+i“' mJ k-i0k^i) = “(f)k Q k 

By using 

9k+i - t? (9k + 0k+8 ) ' 
we have a scheme which conserves 

f S 1 ir 0 te Acr k dS and J £’ n0^Aa k dS 


( 111 . 20 ) 


(111.21) 


under an adiabatic process. Here [ dS is the area integral over the entire globe. 

As in (III. 6), (111.20) is equivalent to 

ir(^- + Wk*v)0k + ^~ [™*k+i(0k +i _0 k) +TI ^k-i(0k"0k-i)] = £^(f ^^k* ^*** 22 ^ 

We define 0 k by 

0k=Tk/Pk X , ("I -23) 
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where p k = p T + mr k - (III 32 ) x p k X gives 

^ (ff + W|c * v )T k - IT 0 k (~- + \V k * v) p k X 

^..cr.-p* e k _!)] = ~Qk / 


from which we obtain 


9 1 r . a i « -| 

ir (0f + W lc *V )c p T k + y^~ |iTO' k+1 Cp(T k+1 “T k ) +mT k _ 1 c p (T k -Tjf-jJj 


= + w«-v)-r 

, i r • 


| — A ^ A ^ ^ A A — 

+ Sa~ L^^+l C P^ r i«+l"P'<: 0k+l) +w *k-l C p(Pk 0k-l _T k-l) 


+ Tf Q v • 


(111.24) 


The left hand side of (111.24) may be written, in the flux form, as 


~(TTCpT k ) + v-(ir\V k CpT k ) + ^ t -(w k+x T k+1 -iw k - l T le _ 1 ) . (111.25) 


8. Total energy conservation and the hydrostatic equation. 


In order that the total energy be conserved under an adiabatic, 
frictionless process, the right-hand side of (111.24), except for TrQ k , must agree 
with ir(u)a) k , where (aja) k is defined by (111.19). First we must require 


(aa) k = a k 


Because (cra) k is defined by (111.16), we have 


j A A R T 

^k _ y\n~~ ^k+1 +1 “ ^k-i CT k-i ) _ lrcr k 
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This is a form of the hydrostatic equation which corresponds to (1.10). We must also 
require that 


•» ^ « A 

c P^"k+x“Pk C P 9k -fx - ^*k “ ^*k +i / 

(111.28) 

and 


Pk c p0k-i“ c P ^k -i ~ ^k-i “®k * 

(111.28) 

Rearranging the terms. 


A A A 

( c P^k+i + ^*k+i) ” ( c P^k +< ^k) ~ Pk c p(Qk+i“0k) • 

(111.29) 

and 


(c P T k +<l> k ) - (cpTk^ + ifc-!) = p* c p (0 k -0 k _ 1 ) . 

(111.29) 


(111.29) and (III. 29)' are analogues of another form of the hydrostatic equation, 
(1.12). 0 k+1 (and therefore 0 k-1 ) is given by (III. 21). 

Replacing k, in (III. 29)' , by k+2 and adding it to (111.29), we obtain 

( C pT)c + 3 + k + 3 ) ” ( C pTic + < £ic) “ C p (Pk+2 + Pk ) (0k+S - ®k) I 
or 

= - e,<p„ x .,- P*> . (111.30) 

(111.30) is an analogue of (1. 11). (111.30) shows that <l? k (at odd levels) is related 

to the neighboring levels as if 0 is linear in p between the levels. 

(111.30) i s usec j f or computing $ at the odd levels. To do so, we need 
to know $ at a single (odd) level, say, k = K. We use (111.27) for this purpose. 

From (111.27), we obtain 

K K DT 

E'^Acr,-^ = E ' Trcr* Act* , (111.31) 

k = i k = i Pk 

K 

which is an analogue to (II. 5). We can write S' <l> k Aa k as 

k=i 
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E'Mok 


k =1 


= V $*( 0 ^- 0 ^) 

k = l 


K-3 


= $K + S' CT k + 1 (4» k - $ k+3 ) 


k=i 


Using (111.30) in (111.32) gives 


K-2 


L' $ k A 0 = + S' cr k+1 c P (p k + 2 - p*) |(9 k+3 + Qj 


k =1 


k=i 


=€> k + s' %r^ic +1 (p^ 2 -p k x ) + <7jc- x (Pk x -Pk- 3 )i9k 

k =i z L 

+ ^*k-:l(Pk X -Pk-3)0k 


+ ¥-n( , -< E K‘) X ) T . • 


Pk 

Substituting (111.33) into (III. 31), we obtain 

R 


$k = $s + S' (~iTCT k — Aa k -c P (o k+1 p k +cr k _ 1 a k ) 
k =1 L Pk J 


T* / 


where 


Pk 


a k = 


|(( E ^ ts -) X -0 for K-2, 

0 for k = K, 

0 or any value for k = 1, 

|fl-(EiLz2) >lN ) fork >3. 

v Pk 


(111.32) 


t\ i T 

jJ T k 

(111.33) 


(111.34) 


(111.35) 
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9. Summary of sections 5-8 . 

We have now constructed a vertical difference scheme which maintains 
the property of the vertically integrated horizontal pressure gradient force, total 
energy conservation under adiabatic and frictionless processes, and conservations 
of 0 and 0 2 , integrated over the entire mass, under adiabatic processes. 


Pressure gradient force (per unit mass) : 
From (111.17) and (111.26), 

RT, 




v$v + o\ 


Pk 


Vir 


(111.36) 


where p k - p T + ira k 


The hydrostatic equation : 

( 111 . 34 ) 

^**k “ c p(tf k + 1 Pk +CT k-i a k)l ^ k . (Ill -37 

i._. L Pv -l 


*=1 


(111.35) 


Pk - 


f i^-') 


a k = 


( 0 or any value 

And from (111.30) and (111. 38) 


for k £ K - 2, 


for k = K, 


for k = 1, 


for k s 3 . 


(111.38) 


^ > k“^ > k+3 c p ( a k+3 ^k+2 Pk^k ) 


(111.39) 
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The first law of thermodynamics 

Using (111.25) for the left hand side of (111.24), rearranging, the terms, 
and dividing by c p , 

~(irT k ) + v-(ir\V k T k ) + [™i k+1 (p k X 0k +1 )-TO k _ 1 (p k x 0 k _ 1 )] 


= * + W lc *v)iT + trQ k /c P , 

which corresponds to (I.26)/c p . 

From (III. 21), 

Pk X 9k+i = i (Tfc + (^^) K Tk+3 ) ( 

P," e'k-! = i ((^ a ) x Tk. s +t, ) _ 


(HI. 40) 


(111.41) 


Among the results given in Secs. 5-8, (111.36-41) are all that we need 
for the main computations. However, for computing moist convection, we also 
need (111.29) and (111.29)', which, with (III. 21), can be written as 


and 


( c pT k+1 +ik +1 ) “ (c P Tk + ^k) = ^ Pk X (0k +3 ”0k) / (111.42) 

( c pT k + $k) - (cpTk-i+ik-^ = •^ E p k X (9k-0k-8) • (111.42)' 


For diagnostic analyses of the results of the computations, such as the 

A A 

vertical transfers of energy and momentum, we may need to know <S> k+1 and T k+1 separately. 

A ^ A 

§ k+1 is obtained by summing (111.27) from k = .1 to k. T k+1 is then obtained from (111.42). 
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IV. VERTICAL DIFFERENCING, PART B 

1. The flux form for the water vapor equation . 

The flux form (III. 5) applied to the water vapor equation is 
J-(itqJ+V-(ir\V k qJ+ T^(TO k+1 q lc+1 -TO k _ 1 q !c _ 1 )=Tr(-C +E) . (IV. 1) 

V/l LMJ J( 

This form guarantees the conservation of total water vapor when there are no 

water vapor sources and sinks. (IV. 1) is equivalent to 

3 1 r “i 

(jjj- + W k *v)q k +^jTO k+1 (q k+1 -q k ) + TW7 k _ 1 (q k “q^ -j.)j 

= (-C+E) / (IV. 2) 

k 

but we must choose q at even levels properly. 

2. Moist adiabatic process — continuous case . 

Consider, first, the moist adiabatic process in the continuous atmosphere. 
Let the air be saturated and remain saturated, and let there be no heating other 
than the heat of condensation.. 

Let q*(T,p) be the saturation mixing ratio. Then, the water vapor 
equation is 

, (iv. 3) 

when condensation is occurring. The first law of thermodynamics is 

c p T = oaa +LC , (IV. 4) 

where L is the heat of condensation per unit mass. 


- 27 - 



IV— 2 


From (IV. 3), we obtain 


, ?al\ dT + - r 


From (IV. 4) and (IV. 5) we have 


C 


1 + 


c P v 9T p 

Substituting (IV. 6) into (IV. 4), we obtain 




where 


dT _ ,8T N 
dT “ “ ( 5p'. ' 


_el _ JL(9a_) 

f — ) = C p C P d P T 

*• '*,7®, 


* ' 


or 


- -[<§?>„- If ! • 


Here 9/3p without the suffix is the derivative under constant horizontal 
coordinates and constant time. 

The corresponding equation with the a-coordinate can be readily 
obtained by using the following relations in (IV. 9): 


(^■ + w *v) p = (^- + \Vv) a -^(~- + \V-v) 1 r| 7 , (from 1.15) 


and 

Then, 


uu = a (g^- + \V*v)ir +TO' . 


<R- + *-V = (i)/(R- + ^ + -[«5?>.-5F 


do 
(1.5) 

3Tx 9T 


(IV. 5) 
(IV. 6) 

(IV. 7) 
(IV. 8) 
(IV. 9) 


(IV. 10) 
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where 


± = 11 . 

3p tt 3cr 


( 1 . 8 ) 


Using the relation 


9q* _ /8q \ . /§g_\ 1L 
^jT - ( df\ + ( #°p 9p ' 


(IV. 8) and (IV. 11) give 

r 3T x 3T _ _ 1 

1 4 

c p v 3 T P 

1 


•wi \ wi _ * r q 3T L 9q i 

^ 9 p n ~ 3 p i i L (dq*) ~'-C P 9 p c p 3 p J 


r _x90 _ L_ 9g^ 

L ^ 3 p c p 3 p J ’ 


1 *£*F> 


(IV. 11) 


(IV. 12) 


From (IV. 10) and (IV .12), we obtain 

r 9T \ ,3 


X 0e + j^3£ 


(i+W.vJJ - o(|f + W-v) ff ir - r ^ a -E- w 

Cp '3T ; 


Equation (1.12) is 


=> p x 0^< c ' T+ *> • 


Defining h and h by 

h =c p T + 3>+ Lq , 

we obtain 


h*= c p T +$ + Lq* , 


x 99 + L 3 g* _ J_ 9 h* 
^ 3 p c p 3 p c p 9 p 


Finally, (IV. 13) may be written as 


(|: +\v.v) a T = (|^W (|: + wv)^ - 


ah’ 

IfL 


(IV. 13) 


(IV. 14) 


: p +L (^ r_) [ 


■* — to 


(IV. 15) 
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3. Moist adiabatic process — discrete case . 

4, 

Let q = q (T k , p k ). When level k is saturated and remains saturated, 

(IV. 2) may be rewritten as 

(|f + W k *v)qt + ^— [Trd k+ i(q Vc+ i-q* k ) + TO k _ 1 (qt-q k _ 1 )j = - c , (IV. 16) 

when condensation is occurring. 

The first law of thermodynamics is, from (111.24), 

(|r + w k -v)T k + ^^r^w-iCp^Qk+rlkH^k-^Tic-pjQk-i)] 


= J“ a k^k(|r + \VvK + ~C , (iv. 17) 


where 


R_L, 

Pk 


From (IV. 16), we obtain 
* 

L _l_, v- 

'9T ' V v 8t 


( i?^ f k <F + ' v « ' 7 )T - + ^ ( l + w ‘ ' v) " 


1 


irAa 


- r^k +1 (qk+i"qt ) + ^k-i(qt-qk-i)]- - c , ( |v - 18 ) 

k u J 


where 


(IV. 17) and (IV. 18) give 


<f£> = <f£) , (f£) M 8 £^> • 

Pk “n Pk “P Tk °Pk Tk 


C =- 




Pk 


- 4t ^ + ™k- 1 (0k-0k- 1 )} 

Pk k 

+ sri^k+iCqk+rqt ) + ^k-^qt-qk-i)}] • (iv.19) 
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Define (3T/3p) nk by 


(f£) 

3 P Bk 


C P C P 9 P TV 

Cp 01 Pk 


Substituting (IV. 19) into (IV. 17), we have 




1 

ttAct, 


[™k+i (Pk X 9 


k+i + c c lk+i“Pk Qk"^ 


* \ 
9k) 


+ ™*k-i(Pk X 9k + —qk 


ic v-. L \ ~] 
- Pv 9*_1- — q k _i) i 


(IV. 20) is an analogue of (IV. 13). 

The coefficient of Tta k+1 in (IV. 20) is 

Pk X (0k + r0k) + 7-(qic + rq?) 

c p 

- Pk ilK^k+a” ®y) + ZT~ fak+i - R k ) 

c p 

= 7-r(c P T k+1 +$ k+1 +Lq k+1 )-(c P T lc +$ lc +Lqt) 

Cp L 

= — (h k+1 - h k ) . 
c p 


Here (111.21) and (111.42) were used. Similarly, the coefficient of Tro k _ 1 
(IV. 20) is 

f<h,*-h M ) . 

C P 


, (IV. 20) 


(IV. 21) 

in 

(IV. 21)' 
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Then, from (IV. 20) we obtain 

(± 

v at 


(lr + \V v ) T k = ^) B!c a k(|r + w k • v)ir 


1 


Cp + L(^-) 




Pk 




(IV. 22) 


(IV. 22) is an analogue of (IV. 15). 

Let us suppose that m7 k+1 is negative and h k+1 >h k . From (IV. 21) 

rW<i*)>p. x (9.-«..i) • ( |v - 23 > 

c p 

From (111.21), 0 k -0 k+1 = §(0k"8k +2 )/ anc * this Is maintained positive or zero 

(dry adiabatically stable or neutral) by the dry convective adjustment which will 

be described later. Therefore, q k+1 >q k when h k+1 >h k . The coefficient of 

•k 

mr k+ in the bracket of (IV. 19) is (?jr- ) p k X (9 k "0 k+1 ) + (qk+a - ^* )/ except 

01 Pk 

fora constant factor. This coefficient is positive. Consequently, the negative 
TrcT k+1 makes a positive contribution to the condensation. (This is not true when 
q k+1 < q k . In that case, the negative ira pumps drier air up from below.) From 
(IV. 22) we see that the negative TTt7 k+1 has a warming effect for h k+1 >h k , 
which leads to a moist convective instability. This instability may occur even 
when no conditional instability exists between the odd levels, which carry the 
temperatures and the mixing ratios. Then the instability is a result of a poor 
choice of q at the even levels and is a kind of computational instability. 
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We may call this "conditional instability of a computational kind (CICK)". 

Similarly, when iw !c _ 1 is negative, h lc _ 1 2: h k is required for 
stability. Therefore, when condensation is taking place at the level k, we 
must choose q k+1 and q k _ 1 which give 


K + ! < h k 
hk-i sh k * 


when *rnr k+1 < 0 


when to k-1 < ® . 


(IV. 24) 


There are three situations for a particular even level which has 


negative to as shown below. 

condensation k 

(a) k+ 1 

no condensation k+2 


h * +1 ^ h k 


(b) 


no condensation k 

— .... 1 ,;;+ ] 

condensation k+2 


bfc+i s b k+ 3 


(IV. 25) 


(c) 


condensation 


condensation 


k 

k+1 

k+2 


bic+3 ^ b k+1 ^ h k __ 


4. Choice of q at the even levels . 

The difficulties of vertical differencing of the water vapor equation are 
not limited to the saturated case. 

For the potential temperature, the arithmetic average, 9 k+1 -^(9 k + 0 k+3 ), 
was used. The integral constraint on 0 2 was maintained and, together with the 
conservation of the average 9, this resulted in a conservation of the variance of 0. 
But the arithmetic average, q k+1 = i= (qk + 9k+a)/ ' s no * a comparable good choice. 
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because, unlike 0, the variance of q is so great that conservations of only its 
lower moments are not effective constraints on its frequency distribution. 

Moreover, the arithmetic average does not guarantee that q remains 
positive or zero . For example, if q k = 0, q k+1 > 0 and in7 k+1 >0, then the 
downward current removes a positive amount from zero . 


The upstream scheme 

9k+i ~ 9k 

9k+l ~ 9k+2 


for 7K7 k+1 > 0 
for < 0 


does not produce this difficulty. But the upstream scheme tends to make q 
homogeneous in the vertical, and thereby produces an excessive condensation 
in the upper levels. 

At present (December 1971) we are testing various choices of q at the 
even levels. 
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V. INTRODUCTION TO HORIZONTAL DIFFERENCING 

1. Distribution of variables over the grid points . 

We now consider the horizontal grid and the way of distributing the 
variables over the grid points. 

Our governing equations are the primitive equations. Under normal 
conditions in the atmosphere (low Rossby and Froude numbers)/ these equations 
govern two well-separable types of motions. One type is low-frequency, quasi- 
geostrophic motion; the other is high-frequency gravity-inertia waves. It is 
known that the energy of locally excited gravity-inertia waves disperses away 
into a wider space, leaving the slowly changing quasi -geostrophic motion behind. 
This process is called "geostrophic adjustment". 

Consequently, there are two main computational problems in simulating 
large-scale motions with the primitive equations. One computational problem is 
the proper simulation of the geostrophic adjustment. The other computational 
problem is the simulation of the slowly changing quasi -geostrophic (and, therefore, 
quasi -nondivergent) motion after it has been established by geostrophic adjustment. 

Winninghoff and Arakawa examined the geostrophic adjustment process 
with various finite difference schemes and found that it depends on how the 
variables are distributed over the grid points. The following discussion is taken 
from their paper. 

★ Winninghoff, F. J. and A. Arakawa, 1972: "Dispersion Properties of Gravity- 
Inertia Waves in Space-Centered Difference Schemes and Their Effect on the 
Geostrophic Adjustment Process. (In Preparation). 
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Dig>ersion properties of gravity-inertia waves in space-centered difference 
schemes and their effect on the geostrophic adjustment process . 

Consider the simplest fluid in which geostrophic adjustment can take 
place: namely, an incompressible, homogeneous, non-viscous, hydrostatic, 

rotating fluid with a flat bottom and a free fop surface. 

The basic equations which govern such a fluid are: 


0.1) 

du - , 3h _ 

dT _fv + s 3: ■ 

0.2) 

dv , , , 3h _ 

dT + fu + 8 37 - 

0-3) 

dh , , /3u , 3v 

dT +h fe + 37 


where t is time, x and y are the horizontal coordinates, u and v are the 
velocity components respectively in the x and y directions, h is the depth of 
the fluid, f is a constant coriolis parameter, and g is gravity. The individual 
time rate of change is 


0.4) 


d _ 3 , 3 , 3 

dt _ 3t U 3x V 3y 


In most of this study we consider the problem with a linearized version 
of these equations. The linearization is done by replacing d/di by 3/3t; and 
replacing h as the factor on (3u/3x +3v/3y) in equation (1 .3) by H, the mean 
value of h. This is justified if the Rossby number is small and the horizontal 
scale is of the order of the radius of deformation. 
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We consider five ways of distributing the dependent variables/ h, u 
and v # in a square grid in space, as illustrated in figure 1 . 



The space finite difference schemes we use for the linearized equations are 
the simplest second-order schemes for each of these five ways of distributing 
the variables. They are: 
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scheme A , 

(1.5) 

|^ - fv + g(6 x h) = 0 , 

(1.6) 

|f + fu + g(6 y h) = 0 , 

(1.7) 

+hF(6 x u) + (6 y v) = 0 

scheme B, 
(1.8) 

|f - fv + g(6 x h) = 0 , 

(1.9) 

|f + fu + g(6 y h) = 0 , 

(1.10) 

|- +H [(6 x u) + (6 y v) = 0 ; 

scheme C , 

(l.n) 

|f - fv Xy +g(6 x h) = 0 , 

(1.12) 

|f + fu xy + g(6 y h) = 0 , 

(1.13) 

|f + H [(6 X u) + (6 y v) ] = 0 

scheme D f 
(1.14) 

|f -*fv xy + g(M>) X V = 0 , 

(1.15) 

|f + FCT X y + g(6 y h) = 0 ; 

(1.16) 

|f + n[(6 x u) + (6 y v) =( 
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scheme E, 

(1.17) 

(1.18) 

0.19) 

where we define/ 
( 1 . 20 ) 

( 1 . 21 ) 


|— - fv + g(6 x h) = 0 / 

|~ + fu + g(6 y h) = 0 / 

fj: + H[(6 x u) +(6 y v)] = 0 ; 

( ^ii * d*"(°i+^ i‘ a i^i) ' 

(a*).. = i (a.., . + a. x .) , 

' M 2 V i+t/| »-£,|/ 


where i and { are the indices of the grid points in the x and y directions. For 

the schemes A through D, d is the grid size d, as shown in figure 1 . For 

scheme E, d equals /T times d, as shown in the figure. This choice of d for 

scheme E will make scheme E have the same number of grid points as the other 

— y 

schemes in a given two-dimensional domain. (6 y a).. and (a ).. are defined 

in a similar manner, but with respect to the y direction. Finally, 

— y 

(1.22) (° = (° *)jj * 

In this study, all analyses with the linearized equations leave the time- 
change terms in differential form. The reason for doing this is that unless an 
implicit scheme is used, we must choose a time interval short enough to satisfy 
the Courant-Friedrichs-Lewy type condition for linear computational stability 
of the wave which has the largest possible phase speed. For the primitive 
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equations of atmospheric motion, this is the Lamb wave. In that case, the 
time interval is adequately small for all other waves, including internal gravity 
waves, and we can then ignore the time truncation error in the first approximation. 
We consider, first, one -dimensional linear equations. The equations 


we use are: 


( 2 . 1 ) 


( 2 . 2 ) 


(2.3) 


3u , , 3h _ n 

ar-fv+9^-0 , 


fr + fu =0 - 


tr +H lr 0 


§fu f H 9fu =0 

0t 2 * ° 9 n 9x 2 


Eliminating v and h, we obtain, 

(2.4) 

I vt) 

If we assume that the solution has a form proportional toe , then 

the frequency v is given by 

(2.5) iyf= l+gH(£) 2 , 


where k is the wave number in the x direction. 

Next, we examine the effect of the space truncation error on the 
frequency. In this one -dimensional case, the space distributions of the 
dependent variables, for the schemes A through E, are shown in figure 2. 
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(A) 

u,v,h u,v,h u,v,h 

• • .. . . i i > 

/-/ i i-H 


(0 

v,h u v,h u v,h 

* — — ■ • • 

/-/ i i-bl 

■(E) 

u,v,h u,v,h u,v,h 

• • 

• / * • t / 

t-j l l+j 

/T 


(B) 

h u,v h u,v h 

• ■ • — • — • 

i-/ i i-H 



(D) 
v u,h 



v u,h 
in 



The forms of the space difference schemes are: 
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scheme A, 



(2.6) 

9u 

9t 

- fv + g(6 x h) - 0 

(2.7) 

9v 

9t 

+ fu = 0 , 

(2.8) 

9h 

9f 

+ H(6 X v) = 0 ; 

scheme 6, 



(2.9) 

9u 

9t 

- fv + g(6 x h) = 0 , 

(2.10) 

9v 

9t 

+ 

c 

II 

o 

** 

(2.11) 

9h 

9t 

+ H(6„u) = 0 ; 

scheme C , 



(2.12) 

9u 

9t 

- fv x + g(6 x h) = 0 

(2.13) 

9v 

9t 

+ fu * = 0 , 

(2.14) 

9h 

9t 

+ H(6 X u) = 0 ; 

scheme D, 
(2.15) 

9u 

dt 

- fv + g(6 x h) = 0 

(2.16) 

9v 

dt 

*«. 

o 

n 

K 

+ 

(2.17) 

9h 

9t 

+ H(6 x u) X = 0 ; 
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scheme E, 

(2.18) |“ - fv + g(6 s h) = 0 . 

(2.19) §f + fu = 0 . 

(2.20) |j: + H(6 x u) = 0 . 

In this one-dimensional case, scheme E is equivalent to scheme A, but with 
a smaller grid size. 

For each of the schemes, we obtain the following frequencies 


(2.21) 

A: 

iff- 

, + ( 

Xf 

.6) 

■ ? 

sin” 2 

kd , 

(2.22) 

B: 

<r> 8 = 

1 +4 

(X\ 

\dy 

2 

sin : 

! (t) - 

(2.23) 

C: 

iff- 

cos 2 ( 

'kd 

v 2 

) + 4 


(2.24) 

D: 


cos 2 ( 

' kd> 

V 2 / 

> + ( 

j sin 3 (kd) , 

(2.25) 

E: 

( r )S= 

1 +2 

(\\ 

\dj 

2 

sin J 



/gFT is the speed of the gravity wave, which is the theoretical maximum 
group velocity of the gravity-inertia waves given by equation (2.5). The radius 
of deformation, X, is defined by^gH/f. The non-dimensional frequency, v/f, 
depends on two parameters, kd and X/d. 

From these frequencies for the gravity-inertia waves in each scheme, 
we can see their dispersion properties. From (2.5) for the continuous case, we see 
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that the frequency of gravity-inertia waves is a monotonically increasing 
function of the wave number k, unless the radius of deformation, X = /gH/f, 
is zero. Then the group velocity, 3v/9k, is not zero unless X= 0. This non- 
zero group velocity is very important for the geostrophic adjustment process. 

The wave length of the shortest resolvable wave is 2d. The corres- 
ponding wave number, k # ' s u/d. Therefore, in examining equations 
(2 .21 )— (2.25), it is sufficient to consider the range 0 < kd < ir. 

Scheme A: The frequency reaches its maximum at kd = ir/2. This 
means that the group velocity at kd = ir/2 is zero. When gravity-inertia waves 
at about this wave number are excited somewhere in the domain (by non-linearity, 
heating, etc.), the wave energy stays there. Moreover, a wave with kd = it 
behaves like a pure inertia oscillation. 

Scheme B: It produces a monotonically increasing frequency for non- 


zero X, in the range 0 < kd <ir. 

Scheme C: The frequency monotonically increases for ^ >% and 
monotonically decreases for Vd < h ♦ For ^ =^, the group velocity is zero 


for all k. 


A. 43 

Scheme D: The frequency reaches a maximum at (g) coskd = \ . More- 
over, kd = ir is a stationary wave . 

TT 

Scheme E: The frequency reaches a maximum at kd = ' 

These results for the one-dimensional case show that scheme B is the most 


satisfactory. However, when X/d is sufficiently larger than 1/2, scheme C is as 
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good as scheme B. For internal gravity waves, it is the reduced gravity which 
matters. When scheme C with a coarse grid is used for an atmosphere which has 
a relatively weak stratification, X/d may not be sufficiently larger than 1/2 and, 
therefore, geostrophic adjustment is poorly simulated by scheme C. This comes 
from the fact that the averaging of the coriolis force, in scheme C, makes the 
shortest resolvable motion behave as if it were on a non-rotating frame. 

Figure 3 shows the case for b = 2 (for which scheme C is also good). 


kd/7r 

F; ? . 3 
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Cahn gave the solution of an initial value problem for which (2.1) to (2.3) 

are the governing equations. At the initial time, he let h = constant, u = 0 and 

v = V in the domain from x = -a to x = a, and v = 0 outside of this domain, 
o 

We obtain the solution, u(x,t), for this same initial condition, by writing 
the solution in the Fourier integral form: 

00 

(2.26) u(x,t)= 2 ~Re J* e'^ X u (k,t)dk , 

— 00 

00 

(2.27) u (k,t) = J e u(x,t)dx , 

— co 

where k is the wave number in the x direction, and u (k,t) satisfies 

(2.28) +(fa +k a g H)u*(k,t) = 0 . 

Equation (2.28) has the general solution, 

(2.29) u (k,t) = A(k) cos(vt) +B(k) sin(vt) , 
where , 

(2.30) v 3 = P(l + \ 3 k 2 ) . 


To determine A(k) and B(k), we apply equations (2.29) and (2.27) at t = 0. 
Then, 


(2.31) 


-ikx 


A(k) - u (k,o) - J e u(x,o)dx = 0 


Cahn, A., "An Investigation of the Free Oscillations of a Simple Current System", 
Journal of Meteorology, Vol. 2, No. 2, June 1945, pp. 113-119. 
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Moreover, equations (2.29) and (2.27) give 

(2.32) ~ 9i^^ = v [~A(k) sin(vt) +B(k) cos(vt)J , 


and 


(2.33) 


9u*(k,t) _ p e -ikx 9u(x,t) ^ 

3t J 9t 


Applying equations (2.32) and (2.33) at t =0, we obtain 

(2.34) B (k) = 1 -i J- e ' ibt (2HM) dx 


t =o 


t=o 


From the initial conditions and equation (2.1), we have 
(2.35) 


^ 9u(x,t) j _ fV for | x | ^a 


t=o 


0 for I x I >a . 


Therefore, from equation (2.34), 


-ikx 


(2.36) 


, a fV e 

B(k) = - r e" ,kx fV dx = - — 

v J o v ik 

-a 


x =a 


x--a 


2fV 


kv 


0 sin(ak) . 


Finally, we obtain 

faV 

(2.37) 

u(x,t) F 

IT 

or 

2fa V 

(2.38) 

u(x,t) 

IT 
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We can also write down the equations corresponding to equation (2.38), 
for each of the finite difference schemes, A through E. These are identical to 
equation (2.38), except for the frequency v, which is given by the expressions 
(2.21) to (2.25) for each scheme. 

We evaluated the integral (2.38) numerically for the differential case 

and for each of the finite difference schemes. We used Simpson's rule, with 

600 intervals in k from 0 to ir/a, to compute u at x = 0 for various values of t. 

Then we obtained h, at x = a, from the equation of continuity. These solutions, 

for a constant x, were calculated for values of t up to 40 hours, at 15 minute 
-4 -1 . 

intervals with f = 10 sec . Since equation (2.38) also enables us to evaluate 
the solution for constant t over a range of x, this was also done for each of the cases. 

Some results of the calculations, with a/d = 1 and x/d = 2, are shown in 
figures 4 and 5. Figure 4 shows the time variation of h for the differential case, which 
approximates the solution obtained by Cahn, and for each of the difference schemes, 
at x = a. Figure 5 gives the space variation of h in the differential case and for each 
of the schemes, at t = 80 hours. 

As we expected, scheme B (together with scheme C in this case where X/d >jg) 
better simulates the geostrophic adjustment than the other schemes. However, even 
scheme B has a difficulty in the two-dimensional case. Figure 6 shows | v | /f as a 
function of kd/ir and 4d/ir, for scheme B, where k and SL are the wave numbers in the 
x and y directions. Again, X/d=2. 
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F '> 5 
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(B) 



It thus appears that all space -centered schemes have some deficiency. 
Because these deficiencies are due to the existence of false low frequencies for 
high wave number motions, we cannot expect that any form of time-differencing 
will avoid these deficiencies as long as we use space-centered schemes at each 
time step. As shown, in Chapter 8, we overcome these deficiencies by using a 
specially designed "time-alternating space -uncentered " difference scheme (TASU). 

2. Two-dimensional nondivergent flow . 

We now must consider the simulation of the slowly changing quasi- 
geostrophic (and, therefore, quasi-nondivergent) motion after it is established 
by the geostrophic adjustment process. 
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We look, first, at the flow which is purely two-dimensional and 
nondivergent. The following is extracted from the review paper by Arakawa* 
(1970). 

The principal computational problems can be illustrated with 
the vorticity equation, 

(13) 9C/3t=J(W), 
where 

(14) J(C,<|0 = @C /3x)(9qj /3y) - (dl£/dy)(d>p/dx) , 

(15) C = v 8 i|) , 
and tji is the streamfunction . 

There are the following integral constraints, among others, on 
the Jacobian: 

(16) Jfctf = 0 , 

(17) P(£7F) = o, 

(18) 9J(C/1>) =0 , 

where the bar denotes the average over the domain, along the boundary 
of which iji is constant. From these integral constraints, we can see that 
the mean vorticity, £, the mean square vorticity, £ 2 , and the mean kinetic 
energy, i(v«|») 2 , are conserved with time. Conservation of these quantities 
during the advection process poses important constraints on the statistical 
properties of two-dimensional incompressible flow. In particular, the 
average wave number, k, defined as 

(19) k 8 = (v 2 q0 2 /(vi|>) 2 , 

is conserved with time, so that no systematic cascade of energy into shorter 
waves can occur. 

If 

Arakawa, A., 1 Numerical Simulation of Large-Scale Atmospheric Motions", 
Numerical Solution of Field Problems in Continuum Physics , (Proceedings of a 
Symposium in Applied Mathematics, Durham, N.C., 1968), Vol. 2, SIAM-AMS 

Proceedings; edited by G. Birkhoff and S. Varga, American Mathematical Society, 
Providence, R.I., 1970, pp. 24-40. 
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If we wish to simulate the statistical properties numerically/ 
we must use a finite difference scheme which approximately conserves 
these quadratic quantities. Avoiding computational instability, in 
the nonlinear sense, is necessary but not sufficient for this purpose. 

Two examples of stable schemes which have a false energy cascade into 
shorter waves will be shown later. 

It should be noted that if we apply Equation (13) to a one dimen- 
sional problem, the nonlinearity will be lost. Therefore, the tests of 
a finite difference scheme for incompressible flow must be made with 
multi-dimensional problems. 

The finite difference approximation for Equation (13) may be 
written, in a relatively general form, as 

(20) Cjk ” £jk ~ AtJj k (£ , ijJ ) , 

where £“ k = (v 2 k <ji) n is a finite difference approximation of £ = V 2 t|J at 
the grid point x = jAx, y = kAy, and t = nAt. Ax and Ay are the grid 
intervals, and At is the time interval. ¥ 2 k and JT jk are finite difference 
approximations for the operators v 2 and J at the grid point x = jAx, y =kAy. 
Hereafter, the suffixes j, k will be omitted unless they are necessary. 

There are a number of schemes corresponding to different choices of £ and 
t|i . For example, £ may be a linear combination of £ n and £ n+1 , such as, 

(21) £* = i(£ n +C n+1 ) , 

•fa 

which is an implicit scheme of the Crank-Nicolson type. Or £ may be a 
provisional value of £, prediction by 

(22) £*=S£ n + aAt/(£ n ,«|» D ) , 

where S and a are 1 , as in the Matsuno scheme, or S is a smoothing operator 

•fa 

and a = as in the two-step Lax-Wendroff scheme. Here J is not 
necessarily the same as J. Or £ may be extrapolated from £ n_1 and £ n , 
as in the second order Adams-Bash forth scheme; that is, 

(23) £* = fc - jC n_1 • 
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The change of mean square vorticity is obtained from (20) , as 

(25) (r +i ) 2 - (rF = At(^ n+i +c )j (c*, «p*y / 

where the bar denotes the average over all grid points in the domain 
considered. Equation (25) can be rewritten as 

(26) (£ n+1 ) 2 - (cW = (C +1 + C n “ 2£*)(£ n+1 -£ n ) + 2At£*J (£*/1>*) . 

To conserve mean square vorticity/ we must properly choose £* and the 
form of J in such a way that the right-hand side of Equation (26) vanishes. 

The first term on the right vanishes if £* is chosen as (£ n+1 + £ n )/2. 

The second term vanishes if the finite difference Jacobian, JT, maintains 
the integral constraint given by (17) on the differential Jacobian, J. 

Similarly, it can be shown that a properly defined kinetic energy is 
conserved if ijj* is chosen as (t|» n+1 + q> n )/2 and J maintains the integral 
constraint given by (18). 

Consider the following three second order, finite difference 
Jacobians: 

A = A x £A y ,|> - Ay£ A x ii> , 

(42) 4 = A y foA x e) - A x («i,A y c) , 

4 = A x (£A y i|>) - A y (£A x t|») , 

where A x /(x) denotes [/(x +d) - /(x-d)]/2d. A y is defined similarly 
with respect to y. It was shown by Arakawa [ 1966] that the Jacobian, 

J, given by 

J=o^ + rJb + P4 / 

(43) 

a + y + P - 1 , 

conserves mean square vorticity if a = |3, and conserves energy if a = y. 

Examples of Jbcobians which have the form of (43) are: 

= + $ 2 ) l 

(44) Js ~ + J 3 ) / 

Jfe = IQk + A) / 

J, = + J 2 + 3U • 

*A. Arakawa, 1966: Computational Design for Long-Term Numerical Inte- 
gration of the Equations of Fluid Motion: Two Dimensional Incompressible Flow. 
Part I. Journ. Computation Physics, Vol. 1, No. 1, pp. 119-143. 
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J 7 Is the Jacobian proposed by Arakawa [1966] . It conserves both 
mean square vorticity and energy. J 2 and J 6 conserve mean square 
vorticity, but not energy. J 3 and J 4 conserve energy, but not mean 
square vorticity. These five schemes are stable. J x does not con- 
serve either quantity, and an analysis similar to the analysis by 
Phillips* [1959], but with the implicit scheme (21), shows that it is 
unstable. J 5 also does not conserve either quantity, but experience 
with numerical tests shows that the instability is very weak, if it 
exists at all. This is not surprising as 2 Js = 3J 7 - J x . Because J 7 is 
a quadratic-conserving scheme the time rates of change of the mean 
quadratic quantities using Js , for given £ and i|», have the opposite 
sign to the time rates of change of the mean quadratic quantities using 

J 7 is the best second order scheme, because of its formal guarantee 
for maintaining the integral constraints on the quadratic quantities. J 7 
is also just as accurate as any other second order scheme. A further 
increase of the accuracy can be obtained by going to higher order 
schemes. The more accurate fourth order scheme, which has the same 
integral constraints as J" 7 , was also given by Arakawa [1966]. 

Numerical tests have been made with the above seven Jacobians. 

In these tests, the initial condition was given by 

(45) q> = '*£' sin(irj/8)(cos(irk/8) + 0.1 cos(irk/4)) , 

and At chosen such that At/d 8 = 0.7. The leapfrog scheme was 

used instead of the implicit scheme. In order to eliminate the gradual 

separation of the solutions at even and odd time steps, which occurs in 

the leapfrog scheme, a two-level scheme was inserted once in every 

240 time steps. The simplest 5-point Laplacian was used. Figures 7 

*N. A. Phillips, 1959: An example of non-linear computational instability. 

The Atmosphere and the Sea in Motion, Rockefeller Press, N.Y., pp. 501-504. 
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and 8 show the mean square vorticity and the energy obtained with 
the seven Jacobians. We observe the expected conservation properties, 
although the implicit scheme was not used. The energy conserving 
schemes, J 3 and JT 4 , show considerable increase of mean square vorticity. 

Figure 9 shows the spectral distribution of kinetic energy, obtained 
by the energy and mean square vorticity conserving scheme I 7 and 
by the energy conserving scheme J 3 , at the end of the calculations. 

The small arrow shows the wave number for sin(irj/8) cos(iTk/8), which 
contained almost all of the energy at the initial time. Although the 
total energy was approximately conserved with 3T 3 , there was a consider- 
able spurious energy cascade into the high wave numbers; whereas with 
F 7 more energy went into a lower wave number than into the higher 
wave numbers, in agreement with the conservation of the average wave 
number, as given by Equation (19). 

Whether the increase of the mean square vorticity is important in 
the simulation of large-scale atmospheric motion will depend on the 
viscosity which is used with the complete equation. A relatively small 
amount of viscosity may be sufficient to keep the mean square vorticity 
quasi-constant in time. However, the viscosity will also remove energy; 
and as a result the average wave number, defined by Equation (19), will 
falsely increase with time. 

On the other hand, the mean square vorticity conserving schemes, 
J s and J 6 , approximately conserve energy, in spite of no formal 
guarantee. This is reasonable, because the mean square vorticity is 
more sensitive to shorter waves, for which the truncation errors are 
large. J 5 approximately conserves both quantities, again in spite of 
no formal guarantees. F 5 , like J x and I 7 , maintains the property of 
the Jacobian, J(£,iji) = - J(<|i,£). T 7 conserves both quantities, with 
only negligible errors arising from the leapfrog scheme. 
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SPECTRAL DISTRIBUTION 
OF KINETIC ENERGY 
AFTER 2,400 STEPS 


♦ 

Jr 


II ^ ^ ^ ^ ^ — — - -- - - - » ■ 

^ (wave number ) 2 



Fig. S 


- 58 - 



V-25 


3* The inertia term in the momentum equation for nondivergent flow. 


Because we are not dealing with purely nondivergent motion, we cannot 
use the scheme discussed in the last section directly as it stands* Moreover, we 
are using the momentum equation and not the vorticity equation. But, the large- 
scale atmospheric motions are indeed quasi-nondivergent. A scheme 
which is inadequate for purely nondivergent motion is almost certainly also inadequate 
for quasi-nondivergent motion. 

Our approach, here, is to first construct a suitable finite difference analog 
to the inertia term in the momentum equation for non divergent flow, and then to 
generalize it to allow for divergence. 

Let us begin with the finite difference nondivergent vorticity equation for 
a square grid in which T 7 of the last section is used. For vorticity, we use the form 


_ 1 


in - (T3 H) * - 4f u ) 


ly tl.l.l-fn *il ^1,1*1 -fu t„ 
d V d d + d d 


(V.l) 


For the grid points shown in the accompanying figure, we define u, v. 


6 y u and 6 x v by 


V.* 5 — 3 — ' 




^1+1, J ” 1*13 


(V. 2) 


(6 > v) 1J (v.3) 
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Then the vorticity given by (V.l) is 

- j«M„ -<«,<-)„) . (v.l)' 

The vorticity equation may be written as 

|j:((6xv) 1J --(6 y u) IJ )=]r iJ (6 x v - 6 y u, <j>) . (V.4) 

Here the symbol J is used for J 7 . Define <|> and q» y by 

(V.5) 

y 

Consider T t j + ^(u, $ ). From a property of the Jacobian, which is maintained 
by U 7 / 

X lf , t |(u/^) = J' I ^ + |(u,‘Ip y +| ud) , (V.6) 

J I/3 - i (u^ y ) = (u # qi T - i ud) . (V.6)' 

— y w y 

Note that (ij» +£ ud) t , + i=(<l> -iud) t ,_i ^ for arbitrary i,j. Then, from 
(V .6) and (V.6)' , 

(6 y J(u,F’))„ - I, iJ ^(u,t')-X lii . i (u,'i y ) 

= r„ (8,u,*) . (V.7) 
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Similarly/ 


(6 X J(v,4i x ) u = J lj (6 x v / il») . 
(V.7)and (V.7) 1 are analogs/ respectively, of 


(V . 7) 1 


37 


J <3y ' ^ ' 


k J <^> 


4 '*’ • 


From (V.7)/ (V.7) 1 and (V.l) 1 / 


(6 x J(v/‘iji x )) 1J - (6 y J(u/il» y )) lj = (6 X v - 6 y u, t|») 

=dJ l3 (C,q») . (V .8) 


We conclude that 


J(u, q» y ) 

for - W • V u 

at the u-points 

and 

J(v, f) 

for - \V’W 

at the v-points 

are consistent with 


J(C,q») 

for - W • v £ 

at the iJj - points 


For purely nondivergent flow, the pressure must satisfy the balance 
equation, which is the divergence equation applied to nondivergent flow. The 
most logioal place to carry pressure is then at the x points in the above figure, where 
the divergence is most simply defined. This configuration is also suitable when we 
are treating pure gravity waves, without (or with small) coriolis force. However, 
we already know that this configuration, which corresponds to scheme C in section 1 , 
is not suitable for the geostrophic adjustment process by the dispersive gravity- inertia 
waves when the radius of deformation is small . 
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For the design of the general circulation model, we are using the 
configuration which corresponds to scheme B. We therefore sacrifice exact 
consistency with the finite difference vorticity equation. The reason that we 
do this is that if the geostrophic adjustment does not operate properly, the 
simulated flow will not be quasi -geostrophic, and, therefore, it will not 
necessarily be quasi-nondivergent. In that case, there would be no point in 
requiring an exact consistency with the nondivergent finite difference vorticity 
equation. 

The integral constraint on mean square vorticity is effective in pre- 
venting a spurious energy cascade, because vorticity is a higher order derivative. 
Because of this, a similar constraint on the inertia term, not necessarily equivalent 
to the mean square vorticity conservation, should (and does) also prevent a spurious 
energy cascade . 

For the differential case, we have 

uJ(u , i|i) = 0 . (V.9) 


In addition, we have 


because 


57 = 0 ■ 


(V.10) 


CV.11) 
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Let iji be carried by those points shown in the accompanying figure. 
Define tp by 


fo) 1+ i, J+ | s i(«P 1J + 4'i +1/ j + \s*i + *i+i,i+i> * 


f 



i.j+l 

U.V 

U.V 

V 

* 

u.v 

U.V 

m 

± 


Here the stream function «j> is defined by 

U l+1,3 + .|. 2d ^1+1 #3+1 + + l ^1+1/3 3 ^ 

v l+^ # 3+l = 2d’^ , l+l / 3+l ^l + l^ ."^1,3 + 1 "^13 ^ * 


(V . 1 2) 


(V . 1 3) 


It can be shown that use of J 7 for J. + x . + i (u, tji) maintains the 

' s 

constraints given by (V.9)and (V.10). Similarly, use of X 7 for J^ +1 , ,+i(v,<|i) 

? 2 

maintains the constraints given by 


v J(v,t|>) - 0 t 

P^ )= °. 


(V.14) 

(V.15) 


The scheme for -\VVu and W-W, which will be given in the next chapter, reduces to this 
Jacobian when the mass flux is nondivergent. 
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VI. HORIZONTAL DIFFERENCING 


1 . The governing equations in orthogonal curvilinear coordi nates. 

Let the orthogonal curvilinear coordinates be £ and 77. The general 
circulation model uses the spherical coordinates/ £ = X (longitude) and 77 = <p 
(latitude). 

Let the actual distance corresponding to d£ be (ds). . Put 


r 4 


and 


(d V • 

(ds) tf n dr > ■ 


(VI. 1) 
(VI. 2) 


Consider a rectangular area element in the £ -17 plane. 





Vw 

- 

/ 


n 




■m 


1 


A£A 77 . 


The actual lengths of the sides are ^ and — . The area is 

m n mn 

Let the component of \V in £ be u and the component of \V in 77 be v. The 


divergence is 




(VI. 3) 


where 6 ^ and 6 ^ are increments in £ and 77 directions/ respectively. In the limit 
as A£/ A 77 -* 0, (VI. 3) becomes 

Similarly/ the vorticity becomes 


r 3 / u, , 3 /V \ 
mn .9£" n + 3tj m j* 


f 3 /Vi 3 /Uv-j 


(VI. 3)' 
(VI. 4) 
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VI-2 


The equation of continuity 


In view of (VI. 3)'/ the equation of continuity (1.16) may be written as 

8 / IT t,9 / U\, 9 / v \ , 9 ,TT<7 




(VI. 5) 


The equation of motion 


where 


The equation of motion (1.22) may be written as 

+ a 0 + (/ + 0 Ik x W + v(i W 2 + <l>) + aavit = Jp , 

S* #CVx\V=mn[|j(^)-|j(^)_. 


The £ -component of (VI . 6 ) is 
9u , • 9u 


9 u . . 9u r . , / 9 / v v 9 ,u»\i 

F + *^-[/ + m n ( 5? (-)-^<-)) J v 

9 .. a . , « . . v . 9 it 


9£ n drj 
+ m — + i v 3+ $) + maa . 


Rearranging the terms# 


9u 


9u 


9u . . 9u 


1 f + mu §f + "v§r^'- /+mn ( V ^n - U ^)] V 


9 l \n 


(VI. 6 ) 
(VI. 7) 


(VI. 8 ) 


+m m +aa 2ii = 


rd<l> , dir i 

_ 9 £ aa 9CJ ' £ 


F t . 


(VI. 9) 


Similarly, 


9v , 9v , 9v , . 9v , r , , / 9 1 9 1 \ 

— + mu r— + nvr- + 5r + / + mn ( v 57 — u-r- - ) 1 u 
9t 94 dr) 9a \_ J \ 94 n dr] m j 


r9$ , dir _ c 

n L^ + " a ^J- F r, 


9ir n 


(VI. 9)' 
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Combining (VI. 5) and (VI. 9), we obtain the flux form for u, 

3 / it \,3 / nu \ . 3 / ttv \ ,3 , tto x 
3t mn U 04 n U 9rj m U 3o mn U 


-r x+ (v 57 i -“|- i )l 

L mn \ o 4 n 077 m /_ 


3 1 Yl 


trv 


. ir f0$ , 0irn_ir c 
+ — XT +aa F . 

n L .04 3^ J mn 4 


Similarly, 


9 / w .A j. 3 , iru x , 0 , irv \ ,3/tro’ \ 

arfe v) + 3 ? ( ~ v) + §5 ' v) + § 7 ( sr v) 


, r / , / 3 1 _ 3 1 

L mn \ 94 n U 077 m 


9 1 \~i 


dr) 


iru 


ir r 9$ . dm it _ 

• - 3 — +aa I = — F 
ml dr) dr) J mn 77 

From this point on, we will consider only those coordinate systems, 
such as the spherical and the cylindrical coordinate systems, in which m and 
n do not depend on 4 • 

From (VI. 10) we obtain the (relative) angular momentum equation 

0 ,71V u. . 9 


3 IT “I 




TTU Ui , d , TTV U x , 0 ,710 U x 

) + sr(— + 


0t'mn m' 04 'n m' 877 'm m 

.r_i Eiif 1t r |_ .£ + oa |_ jli =jl i, 

Lmn mJ L 9 4 mn 9^mnJ mn m 


The first law of thermodynamics 
(1.26) can be written as 

9 / TT T \ . 0 ,7TU T x , 9 /TTV T x , X 9 , TTO -x 

9t ^nn CpT 04^ n CpT) 077 ( m CpT) P 0a ^mn Cp9) 

7T , . U 07T , V 07T ' 


/ 9 / Itv + U dlt , v 07T\ TT n 

mTa ( mr? n 3f mdr)J mn 


(VI. 10) 


(VI. 11) 


(VI. 12) 


(VI. 13) 
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We use the following form: 


3FTu k k * * 

3T +F .*jS,) + 

1 • k+i * It— 1 

+ SS" < s » - s „ ) = 0 - (VI-14) 

Vc 

where 

^=*^2, F. nu^, ShJTo , (VI. 15) 

mn n m 


- 67 - 



VI -5 


and the vertical index k now appears as a superscript. 

For the mass fluxes F and G, the following forms are used: 


*[<“ ■$♦*,>*§ + - ( VL,6) 


where 




4ih 


(VI. 17) 


For the time being, ignore the superior bar operator, which is a linear smoothing 
operator in £ . The form and the role of this operator will be described in 
Chapter VII, section 2. 

Similarly, 

= +1r » ) • (VU8) 


where 




k 4£) 

S ' ' 2 m c 2 


(VI. 19) 


3. The pressure gradient force . 

As in (VI. 10), the pressure gradient force in the ^-direction is 

it r3<f> , 3ir i 

"nl.al + aa 3£j * 

For the first term, we choose the form: 

"nr 3 $ n1c 


(VI. 22) 


" (n 3|) 1+1/j+i "A^A n + w i,j+i^+i,j+i“*i,j+i) 


+ C ir i+V + ,r i,i )( *‘ + V ~®l ,j)j. (V 1 * 23 ) 

Continue, for the time being, to ignore the bar operator. 


- 68 - 



VI-6 


Corresponding to the relation - ^ = - || (ir$) - $ ~] , we 

can rewrite (VI. 23) as follows: 

_ / TT 3$\ k _ _ 1 At? x 

\n9£/ I + £ /J + £ A^Arjn j+| 2 

p 5 " ic ic 5T” “ ’ 

X |_*t+l,3+l *1+1,3 +l“ *1,3 +1*1 ,3+1 i(*l + l,3+l + *l,3+l'(*l+l,3+l~*l,3+l' 

+ *t+l,3*i+i,3“*13 *13-1^1+1,3 + *13> (*1+1,3 "*13 >] * (VL24) 

To be consistent with (VI. 24), we choose the following form for the second term 
in (VI. 22). 

/ ir 3ir \ k _ _ 1 A 77 1 

V n ffa 9£/ A^Atj n , + i 4 

l + £,3+| f 


3tr 


+ (W 1+lf j + Mu) (*1+1,3 “*13 )_ 


where 


*1 + 1,3 +1 “ *1,3+1^ 


,3 " *13>] ' 

(VI. 25) 

k rtl 

7 

v 

(VI. 26) 

p 

r l3 



At each grid point we apply (111.16) to wa in (VI. 25). Adding (VI. 25) 
to (VI. 24), we obtain a form which corresponds to - ~ ( 9 . (?* - "^■(^r) # 

from which we can readily show that the integral properties, discussed in section 2, 
chapter II and section 5, chapter III, are maintained. The momentum generation 
at the grid point i +^, | by the slope of the earth's surface in ^-direction is 
given by 
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/ 


(VI. 27) 


+ (*i + i /3 + 


where 4?^ = (gzs)^ . 

In summary, the pressure gradient force which contributes to 


9t ^ u)l+ i' 3+ £ 


is 


An r k k 

“ n i ^ n i+i /3 +i +ir i / j+i^^ > i+i /3 +i " ^i,j+i ^ 


+ ( ir i+i # J +ir i /3 )( $ i+i /3 - $ i /3 ) 

+ (wI +1/ , tl + W lfJ ti) ( ir i+i /3 +i" ir i , 3 + i) 

+ ((iraa) 1+13 + (wa) 1# ,) (* I+1/3 _ *i,j)] * 

0 k 

Similarly, the pressure gradient force which contributes to g^-(jTv) 1+i J + A 

2 2 

At p k k 

5 " + |j 1T i+i, 3 +i + ir m, 3 )( $ i + 1/3 +i “*i+i, 3 ) 




+ (*j J+1 + * 1 / 3 )(\ 3 + 1 -$i, 3 ) 

+ (Hj H(]+ 1 + Win f3 )( , 1 + Ir3 +i‘ , i«,j) 
+ ( , (Trcra)^ +(wa) x Vir - ir X] . 

' i,i+i i /3 ' ' i /3 +i i /3 /J 


From (VI. 26), 


(rraaf - ir a* 


RT 


ii 


11 13 P 


U 

TT = it — • However ,"JT at the u,v-points has not yet been defined. 


(VI. 28) 
is 


(VI. 29) 



VI -8 


4. Kinetic energy generation ond the first law of thermodynamics . 

The contribution of the pressure gradient force to the kinetic energy 
generation, ^ (JTi u s )* +i/ J+i / is obtained by multiplying (VI. 28) by 
uf.v . + 1 • The contribution of the gradients of and tt at j to the kinetic 
energy generation is 
At; 




X [(* 1 + i , 3 + "i, 3 )(*f+l,J ' + (W 1+l;1 + Wi^) ( ir i+l / i " 1 T i / j)j • 

(VI. 30) 

However, the gradients of and tt inside the brackets of (VI. 30) also contribute 

9 k 

to the kinetic energy generation ^r(J\s^ z ) i + x 3 -i* The combined effect is 

° * 8 ' It 




At** 


*[(*i+i /3 - ^1,J> +((TTora) 1+ i J +(^c.) 1 j )(tt 1+1 J- TT 1W )] . 

(VI. 31) 

As already indicated, the superior bar denotes a linear smoothing operator 
in £. From the form of the operator, we can show that the difference of (VI .31) from 


' 1+ i, 3+ * 


8/ J 8 

k / k 


+ *ijX*i + i,j-*ij) + ((wa)i +1#1 + (iw a )i J )(* 1+lf r*i ,J>] (VL32) 


vanishes when the summation over all i is taken. In other words, 

!R 1 + 1 = 0, 

: T5 


(VI. 33) 
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where 


R, + ^ = (VI. 31) - (VI. 32) , 


(VI. 33)' 


and £ is the summation over all i. Using the definition of F given by (VI. 16), 
i 

we obtain 


(VI. 32) " +i^j “ ^*ij ) 


- l((u “) + (° ~?) ) (M + (naa) ) (tt, + s - ir t , )) 

^ * + £, 3+ | I+ 4, 3 -|r 1+ V 13 

(VI. 34) 

Further, we can show that 
S (VI. 34) = 


- 1 ((u *(i^l ) ((too) +(mro) )(. 1+1 

ov n 1*1 1*1 n i+i .1 -JL/ ' i +1, J 1J ' ' 


1+ £, 3+ £ 1+ ^ 3 ‘ 


" I (f u + < u + ( TOa > )<*» “ w i-i j ^ j 

8V n *-^ 3+ i n »-i, 3 -r v » ‘-v y 7 J 


(VI. 35) 


Similarly, the contribution of the gradients of 4? and it at i to r"*v> 
is given by 

? r<, 


k k k 

j - ',‘*4" G V t 

k 


- i (<"£*■ 
8 


+ (v^-) )(W +(iraa) )(ir + - it ) 

V m i+|,J+£ m i/ M 


S'" 8 


+ ) (W a ) + (™ a ) t j _i ) (■»! , - it, ,_,) j 

8V m i+i, j-i m i-i )-i A ij ' J J 




8/ 8 


(VI. 35) 
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£ (VI. 35) + £ (VI. 35)' , with the equation of continuity (VI. 14), 

gives the contribution of to the kinetic energy generation, 

r ^ TT u 1 .k + i . k— i k 

-[ — + - s .. >]*.. 


(VI. 36) 


Further, following the process which led to (111.18), we obtain the finite difference 
expression for i«a. Using that expression, the thermodynamic energy equation 
(111.40) may be written as 

7 k k k k 

3 k k ^l+i 3 + ^i 3 k ,+Tj-j , 

^n„T u ) + Fj+i , ' 2 ‘ 2 


+ G i j +i ~ J ~ 
l 2 


k k 

1 i J+i 


T, .x. + T t 3 


2 “ G i.3~4 


k k 

T 1 3 +T l 3-1 


/ 2 


l ns fc+1 # k « k+1 k -i-i 

hi J 


I r # Ktl It ,v it-rj. * — 1 * 

+ Ao^ L U ^13^ ®13 “^13 ^13^ 


1 [( \ k LEii. 

C L (lt<Ta ^ 3t 


+ !{(<“ ♦*, . ■ *4 +(u TT^- ‘4, i -j) (H> % ■ + (roa) ") . - *» ) 


+ ((u + ((wa)„ +(wa), _! ,)(!>, 

+ (< v +(v (H» + < TOa C- 1 )K J - V.-i>} 

+ IT Q 1 . ‘ VL37) 

13 13 J 
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5. The water vapor equation . 

The difficulty pointed out in sec. 4, Chap. IV exists even with horizontal 
differencing. However, it is much less serious in the horizontal differencing so that 
we might simply use 

k k 

0 k . Jt +£ hj 




"H, 1 


k k 
+C ?i -1 3 


k J+l +< llJ k + 9ij-i 

+ ~‘~i G - 


Vi 2 


] p.k + i k+i .k-i k-i 

+ Aff7 L 11 4,13 " SlJ qi3 -1 

=7T 1J (- c + e) 


(VI. 38) 


But the arithmetic means shoold be replaced by zero when there is mass flux out 
of the grid point where q is already zero. Possibilities other than (VI .38) are being 
tested (August 1971). 


6. Momentum fluxes. 


The form we chose for gj- (ir ^^2 0 ) + ||(iru ~°u) + At? |— (irv ^-u) 




mn 
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a , Cu) it \ , 1 r 2 r„ k , k , k \ 

aT ( ]Ti^ / 3 + | u ‘ + |/^i ) H 3 l ? * + v + i (Ul+ i' 3+ § Ul+ I, 3+ I 


^ k » 

" ^1 3 +1 \ u i +1 3 +i + U 1 -1 J+i' 
/ J S 8/ 3 3/2 


+ St fl *i <°*-i,<*| +u < k *4,J*i ) ■ "■♦v <u ‘ + i,>*4 +u ‘*l 

1 (u)k+i k+3 k #(u) k-l k , k-S \1 

+ £^?‘ + |, 3+ § (Ul+ i/ 3+ l +Ul+ l, J+ l )_Sl+ l/ J+ l (Ul+ l/ J+ l Ul+ i J+ i J/ 

(VI. 39 ) 



where 


5, +, , +1 - 4 ( F t + 1 J+i + F i + A. ,3 + F i +i_ J + i + + i.^ 


Tl+1 , J+ l 


2 / 


3 / 


3 / 


g [+ im -i( G i+i / 3 + ^ +G i+i / 3 +| ' ~V + i: 
£*, =i(F 1 + i. + F * 

9i j +i “ F i-i 


+ G, , + jl + Gj J+i ) 

/ 3 


■ 1+ 1, m 4v 1 + i /3 + 1 ■ -i+^j+i + ° 1 + V + i + G i + i # J + 4‘ 


_ /3 + 1 - F l+i,J+ 1 + G l /3 +i +G l,J+-l-) 


(VI. 40 ) 
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< u ) 

We have not defined JJ yet. If we put u = const, both in space 


and time, (VI. 39) must become zero. Then we get 





s ‘*4. 


J+l 



) 


. k 


+ fc(3 1+1 m -3 u + g t J+1 -g 1+1 i ) ■ 


.(u)k+i 

rr( s , 

“k 


Acr v 


8 / 


.(u)k- i 

S i+A j+i) 
8 / 8 


0 . 


(VI. 41) 


We derived (VI. 41) rather formally. We can show that (VI. 41) is necessary 
for maintaining the conservation of kinetic energy under a pure advective process. 

From the definitions of 3, g, 3 and g given by (VI. 40), we can show that 

(VI .41) automatically holds if we define 

00 

7T ~ i(T[ 1 + 1 ^+ 1 + rfi+x,j + ITi^+x + TTi^j) • 

(VI. 42) 

= i< S i.V + Vl,i +S l,!*l +s u> (VM3) 


For v, we can use the exactly same form. 


7. Coriolis force. 


See (VI. 10). Coriolis force, plus the metric term which contributes to 


w OT g> is 


[ ! Hz? ' uA « Ar > T n 


(VI. 44) 
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9 

and coriolis force which contributes to (JJ v ) ' s 

_ r r _ u A 3 1 -I 

L mn . ' or? m J 


ITU 


(VI. 45) 


Of course, kinetic energy is not generated by the coriolis force. We must 
maintain the relation (VI. 44) x u . + (VI. 45) x v = 0. We use the following form 
for (VI .44) at the point (I +■§, j +■§). 


]f k k 

8 L^I +1,3+1 + ^i+i / J+i + ^ 1 + 1 , 


k k k 

+ >*>r +7r u > ( C un + C i / 3 


(VI. 46) 


and similar form for (VI. 45). Here 




8 / 2 


2/- S 


(VI. 47) 


When C is constant in i, that is, when the metric term is either negligible 
or constant in i, the zonal average of (VI. 46) is equal to the zonal average of the 
mass flux G, except a constant factor. This is desirable for avoiding spurious 
conversion from Q-momentum to u-momentum. 


- 77 - 



Vll-l 


VII. MODIFICATION OF THE HORIZONTAL DIFFERENCiNG NEAR THE POLES 

1 . Modification of the equations . 

The poles are singular points of the spherical coordinates and the 
velocity components cannot be defined there. Therefore, we let the poles be 
ir-points. 

it at the poles must change as a result of the meridional mass flux, G, 
at all of the points on the latitude circle where the velocity components are 
carried, as shown by the dotted line in the figure. 


7C 



To simplify the computation, we treat each pole as if it were a group of 

points. Each point has index i. At each i we apply the equation of continuity 

3TT 

(VI. 14). After computing and S, at all of the points, we take the average. 

We apply the same procedure to the first law of thermodynamics. It 
follows that those terms in (VI. 37) which are not defined at the poles make no 
contribution . 

In the equation of motion, the momentum flux and coriolis force terms, 
but not the pressure gradient term, must be modified for the points next to the poles. 
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Let the north pole be j = p. We begin by modifying (VI. 42). From 
geometrical considerations/ we let 


(u) x , . 

i+^/P-i ~ Tip + ^vTi +1^-1 + TTi / p-i ' ' 

(VII. 1 ) 

(u) • . • 

(VI 1 . 2 ) 

i + ^/P-^“^p + i^i+i,p-i + ^i / p-i^ • 



{TT^P-i^i/P-i)* |[ 3 { 5 1 + 1 /P .i(u^ /P _i+ ujf + ^ /P _£ 


“ p- l ( u i+i,p- l + u i - i p-i 

' 1 •£' S S' 8 


-95.|,,- a (xf. i ,p- 4 + " 1 . i ,p- 4 )} 


_JL irs (u)!c+1 /^u k+s +u k \ c(u)k-i / k +u k_3 \ 


(VII. 3) 


For energy conservation, we require 
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i 

3 





+ 



(VII. 4) 


3 is not yet defined, but everything else is known. Therefore (VII. 4) determines 
3f , except for a constant part. For the constant part, we choose 


2 r i,P-l 


As for the coriolis force, we let C P = 0. 


(VII. 5) 


2. Averaging the pressure gradient force near the poles. 

For the purposes of illustration, we consider a simple system of equations 
which governs a one-dimensional shallow water wave: 
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3u 


3h 


f + s st = d ■ 


(VI 1.6) 


9h + H — - n 

aF H aF _0 ' 

where the symbols are defined in section 1 , chapter V. Introduce the grid, 
shown in the figure, which is a one-dimensional version of schemes B and C 
(section 1 , chapter V). 

A* il A. u. A- 

— I 1 j 1 1 


(V 11.7) 


I-' i-\ l i+£ t+.i 


d 

Schemes B and C reduce to 



3u +1 

+ J Ov. - h ,> 0 - 

(VII. 8) 


3h t H 

at + d (u <*i' u '-* ) ‘ 0 • 

(VI 1.9) 

Assume that 

A Wi+4) d 

2 

(VII. 10) 


A Ik i d — ___ 

hj =he , where i = /-l . 

(VII. 11) 


Substituting (VII. 10) and (VII. 11) into (VII. 8) and (VII. 9), we obtain 

. kd 

,a sin -=- 

3F + rk (tt) gh = 0 ' 

2 

sin — 

l +rk (-#) HD= 0 • (vl, - ,3> 

2 

A 

Eliminating h, we obtain 
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. kd 


sin 


d 3 u _ , a a / a - 2 \ - 

"3P ” k 0 t Cd ) u ' 


where 


.s = 


gH . 


(VII. 14) 


(VII. 15) 


(VII. 14) is an oscillation equation, and the frequency fi is given by 

. kd 3 

s.n -j 2 

^ = k * c \~dr) • (V,L16) 


For most conditionally stable time difference schemes, the stability 
criterion is given by 

| n | At < e , (VII. 17) 


or 

I c I ^ kd 

-j- sin y< e/2, (VII. 18) 

where e is a constant. For the leapfrog scheme e = 1. sin ^ has maximum 
value for kd = ir, or the wave length L = 2d. To make the scheme stable for ail 
waves, we require that 

|c|At 

-g— < c/2 . (VII. 19) 

Therefore, the smaller the grid size, the smaller must be the value of At. 

Because the meridians converge to the poles, the grid-size in the 
^-direction becomes very much smaller than the average grid size over the 
globe as a whole. As a result, an extremely small At must be used to assure 
stability. 
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We might be able to use different At for the different latitudes, but 
that procedure would be very complicated. The usual approach is to decrease 
the number of grid points on the latitude circle as the latitude approaches the 
pole. However, it will be extremely difficult to design a space difference 
scheme, for such a grid, which will have the integral constraint given by (V.10). 
Another procedure is to keep the regular spherical grid, but use a larger space 
interval to compute the finite difference quotient; but this decouples each grid 
point from its neighboring points and introduces spurious computational modes. 

The method used in our model is to smooth the £ -component of the 
pressure gradient force and the divergence in the ^-direction. Let the smoothing 
operator modify the amplitude of the pressure gradient and the divergence by 
the factor S(k). (VII. 12) and (VII. 13) are replaced by 

do - / s ' n 

5F + ik ( — i«T“) 9 s(k)h = 0 ' (VH.20) 

T 

A . ,kd» 

jl _ sin(-y) 

3t ik (“IT“) HS(k)G = 0 • (VII .21) 

T 


The stability criterion (VII. 18) then takes the form 


c I At 


If we choose S(k) such that 


2 — sin £§r) s ( k ) < f • 


sin(~) S(k) ^ , 


(VI 1.22) 


(VI 1.23) 
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where d is a specified standard length , then 


! G j At 
d 


sin(~) S(k) 


I c l At 

d* ' 


so that the criterion for stability becomes 

| c | At 
d* < 2 


(VI 1.24) 


rC 

This criterion depends only on the specified standard length d . 

In the model/ we let d be equal to the latitudinal grid size d = 

Tj n 


which is a constant. The longitudinal grid size is d^Q) h 


i) = M. 


m, 


Then 


d *(i) /k » 

S(i/k) =f-- / sin(-d Q)) / 

V * 


(VI 1.25) 


when the right hand side of (VII. 25) is less than 1 . Otherwise S^‘,k) = 1 . 

To do this, we expand the zonal pressure gradient and the zonal mass 
flux into Fourier series and reduce the amplitude of each wave component by the 
factor S(j,k). These are the bar operations shown in chapter VI. The smoothing 
in the mass flux given by (VI. 16) is chosen to maintain the energy conservation. 


It should be noted that this smoothing operation does not smooth the fields 
of the variables, because it is simply a generator of multiple point difference 
quotients. For the example that is given above, the solution of (VII. 20) and 
(VII. 21) is a neutral oscillation. 
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Vlll-l 


VIII. TIME DIFFERENCING 

1 . Shortcoming of space -centered schemes . 

As was pointed out in section 2, chapter V, all space -centered schemes 
appear to have some deficiency when the time change is in differential form. 

This situation is not modified by any time differencing scheme, if the space- 
centered differencing is used at each time level. Every space-centered difference 
scheme introduces averaging in either the pressure gradient force or the coriolis 
force, and this is the cause of the trouble. 

To illustrate this difficulty in a most dramatic way, the following 
experiments were done by Winninghoff and Arakawa. They integrated the two- 
dimensional shallow water equation, keeping the coriolis force and the advection 
terms but no physical dissipation terms, and using scheme B and different forms 
of time differencing. They prescribed a point mass source and a point mass sink 
which were 10 grid intervals apart. The figure on the next page shows the height 
of the free surface, h (when the average value of 1 km is substracted) , around 
the sink point, for four schemes of time-differencing. All of the schemes, but 
one (the TASU-Matsuno scheme), show a false checkerboard pattern of the free 
surface height instead of the monotonic slope toward the sink. 

In the case of stratified flow, heat sources and sinks will have the same 
effect on the temperature field as the mass sources and sinks have on the free 
surface height field in these experiments. 
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TASU-Matsuno 


VIII— 2 vm “ 2 



VIII— 3 


Vltl-3 


2. The TASU scheme (time-alternating space -uncentered scheme.) 

Consider the two-dimensional, shallow water equation. For each of 
the directions, x and y, we use a one-sided space difference at one time level. 
But to obtain an overall accuracy comparable to the centered difference, we use 
the one-sided space difference at the opposite side at the next time level. 


The TASU (time'Xilternating space-uncentered scheme): 
At even time levels, 

(|^ij * d 




j- upper-uncentered (VIII. 1) 


0” * 


V 1 **,^* * j+i _ 


At odd time levels, 

^ 1+ §' J + i * d^ 1+1 ,* ~ h V ^ 


l?>« - 


d ( V l-i,3+i‘ V l 

a s' s s' s 


*9yW' J+ ! * d ■ h i,j ? 


| right-uncentered (VIII. 2) 


j- lower-uncentered (VIII. 3) 


j- left-uncentered (VI 1 1.4) 
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V1II-4 VIII-4 

k, k A, 

k 

t+'.j 

If there is no advection and no coriolis force, that is, if we are dealing with 
pure gravity waves, it can be shown that a simple use of the Euler scheme with 
the above time alternation is stable for small At. With advection and coriolis 
force, the simple Euler scheme becomes unstable, so that we have to combine 
the above alternation with some other scheme. 

A scheme which combined the TASU scheme with the Matsuno scheme 
was tested for the example described in section 1 . To explain the procedure, 
we write the equations symbolically in the following form: 

/(A) (VIII. 5) 

at 

The leapfrog scheme is 

= (An) 

2At V A } ' 

The regular Matsuno scheme is 

A(n+i)*_An 

At ' > ' 

An +1 - An . * 

A a . A- = /c( aU+i> ) . 
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VII 1—5 


VI 11-5 


Here / c isa space-centered difference scheme for f . 
The TASU-Matsuno scheme is 


Time Space 


A (n + l) - A n _ 
At 

/ C (A“) 

EULER (Forward) 

centered 

A »+i - A n _ 

At 


Backward 

upper-right 

uncentered 

A( » + a)* - A n +1 

• = / c ( A*«) 

EULER (Foreword) 

centered 

At 

A n + 3 - a d + 1 

= /„ l(A‘* m) *) 

Backward 

lower-left 

At 


uncentered 


The results obtained by these three procedures, as well as by a version of the 
Lax-Wendroff scheme, were shown in the figure on page VI 1 1-2. 


3. Time differencing of the model . 


The upper-flux in the £ -direction is defined by 
(Fl Supper = 


Similarly, 


^i+i^lower + §,J-iK+ V +iri J> ' 

(Gi ,j + iWt = ^ (v ^W' 3 + i (l v +i+1Tij) ' 

k 

^ Gl ,J + ^left = !(v ^)t-i A3+ i K fJ +i +ir 13 ) • 


(VIII. 5) 


(VIII. 6) 
(VIII. 7) 

(VIII. 8) 
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VI II— 6 


Corresponding to (VI. 28) and (VI. 29), we also use the upper, lower, right and 

left pressure gradients. Also, in the right hand side of (VI. 37), only one of 

the (u 4ll) or ( v terms within each bracket is picked up, with g replaced 

4, when we want the uncentered expression. 

4 

We must use uncentered expressions at the same side everywhere in 
the system of equations at a given time, level. Otherwise, an instability may occur. 

The actual time-marching procedure used in the model is primarily 
the leapfrog scheme with a periodic use of the TASU-Matsuno scheme. At the 
present time (December 1971), the TASU-Matsuno scheme is used every fifth time 
step. 
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IX. LARGE-SCALE PRECIPITATION, DRY AND MOIST CONVECTION 
1 . Large-scale precipitation. 

k *k 

Large-scale precipitation occurs when - q is the excess mixing 
ratio; but not all of this excess condenses, because the temperature, and 
therefore q l3 , increases as a result of the condensation. A first guess of the 
condensation is taken as. 


CAt = 


_ k _ _*k 

]+ r^ 

C p O I plc 


(IX. I) 


CAt is subtracted from q k , and LCAt/c p is added to T k . This process may be 

U i 3 

iterated for a better accuracy at the given step. 


2. Dry convective adjustment . 

k k+3 


When 0^ < 0 U for any odd k K -2), we assume that subgrid-scale 

k k+3 

dry convection occurs. We modify T tJ and T J3 in the following way: 


(AT^Ac* = 

- (AT U + ) Acr^s 

(IX. 2) 

k 

k+2 


k + at u 

k+3 AT j3 

(IX. 3) 

/ it \X 

9 >i + / 

O 

0**0 



This temperature change at levels k and k + 2 may cause a new unstable lapse 
rate at a neighboring interval. If it does, then the above adjustment is applied 
to that unstable interval. This process is repeated until all of the intervals are 
stable. In the three level model two steps are sufficient. 
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3 . Parameterization of cumulus convection . 


Consider an ensemble of cumulus clouds. Here, the cloud is defined 
as the saturated portion of the air. Let be the fractional area covered by the 
ith cloud in a horizontal cross-section at level z. The total fractional area 


0 

0 

0 


c 


<Y^ 

O 


<?. 


covered by all clouds is a = Sa., where E is the summation over all clouds 

i i 

in the unit area A. 

/ 

As a basic continuity equation, we let 


V‘(P\V) + (PW) - 0 


(IX. 4) 


where p is a function of z only. 

The total vertical mass flux in the cloud ensemble is given by 

M(z) = E J pw da . (IX. 5) 

i a i 

We define the large-scale vertical mass flux by 

to (z) = J pw da= pw . (IX. 6) 

A 

M-to is the net downward mass flux in the environment. The schematic cloud 
in the accompanying figure represents the ensemble of clouds. 



unnmrnjrrnrrT 
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IX -3 


Define h by 

h = c p T + gz + Lq , (IX. 7) 

where c p T, gz and Lq are, respectively, enthalpy, geopotential energy, and 
latent heat, h is approximately conserved with respect to an air parcel, or 

gp = 0 . (IX. 8) 

Combining (IX. 8) with (IX. 4), we have 

(ph) + V *(p\Vh) +|^(pwh) = 0 . (IX. 9) 

We must consider entrainment and detrainment layers separately. 

Inside the clouds, in an entrainment layer, and immediately outside of the 
clouds, in a detrainment layer, we expect strong turbulent mixing. In the 
entrainment layer, we integrate (IX. 9) over an area slightly larger than the area 
ffj , at the boundary of which turbulent mixing may be ignored. Then we obtain 

= 0 • (IX. 10) 

where 

Mj = pWj da . (IX. 12) 

Here, the vertical transport of h by the internal structure of the cloud is neglected. 
3Mj/9z +p9a./3t is the entrainment rate of environment air into the cloud, which 
may not be steady. 

In a detrainment layer, we integrate (IX. 9) over an area slightly less than 
the area a l at the boundary of which turbulent mixing may be ignored. 


h = — f h dcr , 
1 cr, J cr, 


(IX. 11) 
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Then we obtain 

§f(p h i°i>-(i^ + ptn )h i + |l (M i h i) = 0 • ( ,x - 13 ) 

If we assume that the individual clouds are alike, then the summation 
of (IX. 10) over all clouds gives 

(P h c ff ) " (|f + P |f) h e + §7 (Mh c ) = 0 , (IX. 14) 

where h c = hj. Similarly, from (IX. 13), 

|f (ph c <r) - (— - + p|^)h c + |^(Mh c ) = 0 . (IX. 15) 

Integration of (IX. 9) over the environment gives 

|f (ph,(l-a))+ (^-+ P |f) h 0 - |(Mh s ) = - (pwh 0 ) - v^pWhJ , 

(IX. 16) 

for the entrainment layer. Here \V is the large-scale horizontal velocity. For 
the detrainment layer, 

|f(pb 0 (l-a))+ + p|f0 h c - (M h e ) = -~(pwh 6 ) - V-(pWh e ) . 

(IX. 17) 

We assume that a is much smaller than 1 . Then we obtain the following 
approximate equations J 


for the entrainment layer: 

from (IX. 14), ~ h, + ^ (Mh c ) = 0 , 

from (IX. 16), 3 , u .. 3h* _ 3 , — , \ , - , * 

9f (P h e) - -07( P wh e )- v(p\Vh e ) , 


(IX. 18) 
(IX. 19) 


See the last paragraph of this section. 
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i 


for the detrainment layer: 
from (IX. 15), M 

from (IX. 17), _9_ 

9t 



r 


(ph s ) +(h c -h 9 ) -M 



(IX. 20) 

-^■(pwhj- V.(p\Vh e ) . 

(IX. 21) 


From the definitions, h = h c a + h e (l - a). Because we are assuming 
that a 1 (and h c ~h e ), we approximate h s by h. Then, 


for the entrainment layer: 

'ir 5 ' + s' (Mh = ) = 0 - (lx - 22) 

for the detrainment layer: 

M^=0. (IX. 23) 

In the environment no condensation occurs. Separating h into s and Lq, where 
s = c P T +gz, we obtain 


for the entrainment layer: 


If (ps) = M ~ (pw s ) - v .(p W s ) , 

(IX. 24)** 

|f (pq) =M || - ~(pw q) - V-(p \V q ) , 

(IX. 25) 

for the detrainment layer: 


|f(ps ) = - (s c -s ) |^ + m|| “ |f(pw s ) -v.(p WD , 

(IX. 26)** 

|r(pq) = " (q c -q)^- +M ^ “ ^(pwq ) - V-(pWq) , 

(IX. 27) 


(IX. 22-23) are diagnostic equations which determine h c as a function of height, 

and (IX. 24-27) are prognostic equations for the large-scale temperature and 

*lt is also assumed that there is no evaporation of liquid water in the environment. 
See the last paragraph of this section. 

**More exact equations are obtained by replacing s by the potential temperature, 9 • 
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water vapor fields. 

Condensation in the clouds is given by water vapor inflow from the 
environment, so that the environment water vapor loss, due to the flow induced by the 
clouds, is the amount of condensation. From (IX. 25) and (IX. 27), this loss is 



For the entrainment layer, the convective warming of the environmental 
temperature (which is the large-scale temperature) is through M This is also 
approximately true for the detrainment layer, because s c -s is usually small there. 
The warming effect through M ^ represents the adiabatic warming due to the 
subsidence M in the environment, which compensates the upward mass flux M in 
the clouds. (See the figure at the bottom of page IX -2.) It should be realized 
that the heat of condensation released in the clouds is used for maintaining the 
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excess temperature of the clouds against adiabatic cooling and entrainment of colder 
and drier air, and, therefore, is used for maintaining M. The heat of conden- 
sation is not directly used for warming the environmental air. (However, whether 
net environmental subsidence exists or not depends on the sign of M -!T\.) The factor 
Tj at a certain level in CISK models (for example, Ooyama (1969)), may now be 
interpreted as the ratio of M at that level to the large-scale upward mass flux at the 
top of the boundary layer which is produced by the mass convergence below. 

Our problem is to find M(z). If M(z) is somehow determined, we can find h c 
from (IX. 22) and (IX. 23). Then we can find s c (and, therefore, T c ) and q c from 
(IX. 29) and (IX. 30). The condensation can be calculated from (IX.28). The temper- 
ature and mixing ratio of the large-scale fields are predicted by (IX. 24-27). 

The above discussion indicates that relating the mass flux M(z) to the large- 
scale fields must be, at least logically, the central core of a cumulus parameterization 

scheme . 

It is plausible to assume that cumulus convection adjusts the environment in 
such a way that the energy supply from the environment to the cumulus convection is 
eventually terminated and M(z) becomes identically zero, unless a counteracting 
modification of the environment by large-scale processes exists. Here, radiation and 
sensible and latent heat supply from the earth's surface are included in the large-scale 
processes. When modification of the environment by large-scale processes exists, the 
environment may gain energy in the form available for cumulus convection. We then 

*Ooyama , K., 1969: Numerical simulation of the life cycle of tropical cyclones. 

Journ. Atm. Sci., 26, 1, pp. 3-40. 
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consider a quasi-equilibrium of fhe cumulus ensemble in which the large-scale 
processes, acting as forcing functions, are balanced by the convective adjustment. 

M(z) can be determined when we assume this balance. The feasibility of a parameter- 
ization of cumulus convection crucially depends on the existence and uniqueness of 
such an equilibrium state. When M(z) is found, we can estimate the characteristic 

cloud size in the entrainment layer provided that the relation rjr ~ = — holds, at 

M oz r 

least approximately. 

Based on the considerations outlined above, A. Arakawa (1969)* proposed a 
parameterization of cumulus convection for a three-level model. The current general 
circulation model (December, 1971) uses essentially the same parameterization. 

Recently, Arakawa (1971) presented a new parameterization, which is physically 
more realistic and is applicable to any number of levels. The assumption that all clouds 
are alike is abandoned and, instead, a spectral distribution of cumulus convection is 
considered. Evaporation of detrained liquid water in the environment and interactions 
of the cumulus convection with the subcloud layer are taken into account. We are 
in the process of testing this new parameterization in the general circulation model. 


K 

A. Arakawa, 1969: Parameterization of cumulus convection. Appendix I, Numerical 
Simulation of the General Circulation of the Atmosphere. Proceedings of the WMO/ 
IUGG Symposium on Numerical Weather Prediction, Tokyo 1968, pp. IV-7 to IV-8-12. 


A. Arakawa, 1971: A parameterization of cumulus convection and its application 
to numerical simulation of the tropical general circulation. The Seventh Technical 
Conference on Hurricanes and Tropical Meteorology, December 6-9, 1971 , Barbadoes. 
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X. SURFACE FLUXES AND PREDICTION OF GROUND CONDITIONS 


1. Surface friction. 


** 


The surface velocity, \V S , is estimated by a linear extrapolation of \V, 
with respect to a, from levels 3 and 5 to the surface, which is level 6. The 
surface stress is 

T s = - pc D I \V S I \V S . 

When the surface air temperature, T s , is equal to the ground temperature T g , 
the "surface layer" is neutral, and then 

over open ocean: c 0 = 0.001 x (1 + 0.7 | \V S | ) , 

with c D =0.025 as the upper limit. 

over land, ice or snow: c 0 = 0.002 + 0.006 x 2 ^ /5000, 

where z s is in meters. When the surface layer is not neutral, over all surfaces, 

1 


= (c D ) 


neut. 


l-7.0rpr, for AT h (T K -Ts) < 0 , 


C[j (c D ) 


neut. 
-1 


fAT 


) 


for AT > 0 , 


0+ JWs I s 

where AT is in °C and | \V S | is in m sec 1 . When computing c D , AT is taken 
as the average of the newly computed AT and AT at one time step earlier, in 
order to avoid oscillations in time. 


2. Surface sensible heat flux 


** 


The surface sensible heat flux is 


Fs c p p C D Wsltf, - Ts ) . 


The formulations in this chapter were mainly done by Dr. Akira Katayama. A more 
detailed description will be published as a technical report. Numerical Simulation of 
Weather and Climate, Department of Meteorology, UCLA. 

**See the footnote on the next page. _gg_ 



X-2 


To compute T s / we assume that F s is equal to the flux at the top of the surface 
layer, which is taken as 


" c pP K((T 5 -T s ) - ( T 5 -T s ) crit )/(z 5 -z s ) . 

2-1 

We are currently using K = 10 m sec . We choose (T 5 - T s ) . in the following 

criT • 

way: 


and 


05 - T s ) , - 05 - T s ) dry adiabatic, when r s = 0 , 


05 -T s ) . - 0 b - T s )moist adiabatic, when r s =1 . 

crir 


Otherwise 0 S - T s ) c ^ >s linearly interpolated between 0 S - T s )^ and 

05 - T s ) with respect to r s , where r s is the surface relative humidity, 
m .a . 


3. Evaporation . 

For ocean, snow and ice, the evaporation is 

E s = pc 0 | W s | (q*0 K ) - q s ) . 

q s is determined in a way similar to determining T s , but without (qg - q s ) . 

cr i T 

The evaporation from bare land is taken as 

E s -P E S/ p 

where 

E s,p = P c D |W S [ (q*0 B ) ~ q s ) 

= P c 0 | W s 1 0 g - T s ) + q* - q s _ 

is the potential evapotranspiration, and p is a coefficient which depends on 

the wetness of the ground . 

* * 

We are in the process of modification by introducing an explicit planetary 
boundary layer following the line proposed by J. Deardorff, Mon. Weather 
Rev., j00, (1972), pp. 93-106. 


- 100 - 



X-3 


We let w be the amount of water in the ground, per unit mass of ground, 
which is available for evaporation, where w^ is its maximum possible value, 
and we define the ground wetness by w 1 = w/w^. For relatively large w' , 
capillary forces are sufficient to carry the water to the surface for the potential 
evapotranspiration, so that we let 


0 

x 

0 0.5 1 

P = 1 when w‘ ^ 0.5 , 

P = 2w’ when w 1 < 0.5 . 

However, when q (T, ) <q s , that is, when there is negative evaporation (dew 
deposit), P = 1 . 



4. Prediction of ground conditions. 

The table on the following page shows, by a check mark, those variables 
which are predicted for the different types of ground condition. 


Ground temperature 

For sea ice, we assume a uniform thickness of 3 m, and calculate 
the heat conduction through the ice from the ocean below. For the prescribed 
permanent land ice and for the calculated snow over soil or ice, heat con- 

duction from below is neglected. The heat capacity, C, is given by — , 

0 ) 
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X-5 


where X is the thermal conductivity, c is the heat capacity per unit volume, 
and a) is the angular frequency of the diurnal change. 

The prediction equation for T g is 

C S s - R s F s - LE S + Q h , 

where , at the ground surface, S s is the net downward flux of solar radiation, R s 
is the net upward flux of long wave radiation, F s is the upward flux of sensible 
heat, LE S is the upward flux of latent heat, and Q M is the heating (or cooling) due 
to the freezing of water (or the melting of snow or ice). R s , F s and E s depend 
on T g . . That part of R s which is emitted by the ground depends on T ? and is 
proportional to T ( 4 . 

We solve the prediction equation for T g by the backward implicit method. 

n + i n 

We do this by linearizing the equation with respect to the unknown, AT t = T g -T t , 
where n is the index of the time level- 

Snowfall, snow melting and ice melting 

When T s < 273. 1°K, we let precipitation take place as snowfall 

(or snow showers). The mass of snow over land or ice increases with snowfall 
and decreases by evaporation and melting. Melting occurs when the predicted 
temperature, without Q H , is higher than 273. 1°K. Then T g Is adjusted to 273. 1°K 
by adding negative Q H , with a corresponding reduction in the amount of the 
snow mass. 
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Similarly/ melting of the prescribed land and sea ice occurs, 
with the adjustment of T, to 273. 1°K by adding negative Q M ; but the mass 
of the ice, not being a time-dependent variable in the model, 
never melts completely away. 


Ground wetness P £ 



Fv 


For bare land, 

§H3 = p _ F - F - F 

Qj. ~ c r v r H / 

where 

-2 

m is the total water in the surface soil layer (g cm ) 

-2 -1 

P is the rate of rainfall (g cm sec ) 

-2 -1 

E is the surface evaporation ( g cm sec ) 

F v is the downward drainage of water through the lower boundary of the 

surface soil layer which holds the available water 
F h is the horizontal drainage of the available water in the soil, including 
the runoff at the surface 

Also, w = m/ph, where p is the bulk density of the soil and h is the thickness 
of the layer which holds the available water. 
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Then 




From the definition of the ground wetness. 



9w' _ P.-E -(F h + Fv) 

3t w m ph 

The total drainage runoff, F H + F v , depends on the rate of rainfall and the 
wetness* The functional form currently used is 


F h + F v = (P 3 + D 3 )^ - D , 

where D is the water deficiency in the soil layer defined by 

D = (1 -w')iw ph . 

' m r 


When w' = 1 , all of the rainfall runs off. We are letting % ph, the total 

-2 

maximum available water in the ground, per unit horizontal area, be 10 gr cm . 

The ground wetness is also predicted for the soil under the snow. Both 
rainfall, which penetrates the snow, and snow melting contribute to an 
increase of the available water stored beneath the snow. 

The intersticial water in the soil may freeze, and become intersticial 
ice, when T g <273. 1°K. This type of ice is carried as a prognostic variable. 
When T f becomes larger than 273. 1°K, it melts. This freezing and melting also 
affects the ground temperature. 
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XI. RADIATION* 
1 . Long-wave radiation . 



T<*> 


T 


r 3 


The net flux (upward positive) at the level z is given by 

T 

+ J* T(u*- uf, TjTrgydT + f 't(u*- u*, T)irg|dT , (XI. 1) 

O T * ‘ 

* Z 
where u is the effective amount of radiating substance integrated from the 

earth's surface to the reference level z, t is the mean transmission function, averaged 

over the whole frequency range with the weight dirBy/dt, ttB v is the black body radiation 

at frequency v, and ttB = aT 4 , the black body radiation integrated over the whole 

frequency range. 

A typical vertical profile of t, as a function of crT 4 , is shown in the 
1.0 

accompanying figure, 

T 

f 

T 

° <*T<% t£ aTj 

The area below the curve is the integral (XI. 1). 

The formulations in this chapter were done by Dr. Akira Katayama. A more detailed 
description will be published as a Technical Report, No. 6, Numerical Simulation of 
Weather and Climate, Department of Meteorology, UCLA. 



- 106 - 



XI— 2 


_ _ c* r* 

The variation of t Is mainly through u - u z or u z - u , for 
temperatures higher than about 220°K (= T c ). For T >T C , we can use a more 
or less arbitrary average temperature T for t(u - u z , T) or t(u z - u , T). 

We may approximate (XI. 1) by 

R Z =A+B + C, 

where the sub-areas A, B and C are given by 

i _ i _ j 

A = aTg T(u m - u z , T c ) , 

B = (aT 4 - aTt ) t(u* - u* , T ) , 
crT 4 

C = J t( | u - u z | # T ) d(aT 4 ) . 

aT 4 

00 

We choose T = 260°K.‘ T is the mean transmission function, averaged over the 
whole frequency range with the weight trB v . 

A straightforward computation of the integral C requires very high 
resolution near the reference level because of the sharp maximum of t. 

The computation of flux in our general circulation model is done at 
the even levels (0,2,4, 6). Suppose we are computing the flux at level 4. For 
the layer 0-2, we may use the simple trapezoidal rule. But for the layers 2-4 
and 4-6, we need to take special care. We assume that the average of t 
in the layer 2-4, for example, can be approximated by 0 +m 24 T 2 
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The optimum value for m was determined empirically. We found that we 
can safely use a constant for each layer except for small regions such as high 
mountains and for extreme situations such as very high or very low mixing ratio. 

For the vertical spacing given in the figure, the constant optimum values for m are as 
shown below. ~~~ ™“ 


T 


to 


Co 

VO 


* 

6 


o 

T) 


O 

oo 


too mb 
400 mb 


-ftn to o mb 

— e 0 ' j) *° /ooo mb 


A more general optimum value of m as a function of various parameters such as pressure, 
temperature and mixing ratio is being tested. 

For the transmission functions, the following empirical formulas are currently 


used; 


for H„0: t(u*, 260°K) = 


1 


t(u*, 220° K) = 


1 + 1.75 u 

1 


*0.416 ' 


1 + 3.0 u 


*0.408 ' 


for CC^: The effective amount of CC^ between the levels 

p 1 and p (p 1 > p) is 

u c o (p* *P) = 127 P / poo = 1000 mb 

3 POO 

(unit: cm-NTP) 

T co 2 ( u cog ) — 0.930 ~ 0.066 log^Q u co 2 • 


In the 3-level model, cloudiness is either one or zero. For large- 
scale (non-convective) clouds, and T < - 40°C, we have the greyness factor a = 0.5. 
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2. Solqr radiation 

Following the suggestion by Joseph (1971), solar radiation is divided 

into the 

scattered part, Sq =0.651 S 0 cos C , 

and the 

absorbed part, S® = 0.349 S 0 cos £ , 
where S 0 is the solar constant and £ is the zenith angle of the sun. 

The absorptivity for the absorbed part is 0.271 (u sec • 

We define the function A(X) by 

A(X) =0.271 X 0,303 . 

The albedo of the atmosphere due to Rayleigh scattering is 

p 

Oq = 0. 085 - 0 . 245 log . n cos £ ) , 

I u r o0 

where P 3 is surface pressure in mb and P 0 o is 1000 mb. 

For the albedo and absorptivity of clouds, the values used are those in 
the table below, as given by C . D . Rodgers (1 967). 


Cloud 

Scattered part- 

Absorbed part 


Type 

albedo (Ri 1 ) 

albedo (Ri) 

absorptivity 

(Aj) 

n 

u c t 

1 

0.21 

0.19 

0.04 

0.01 

2 

0.54 

0.46 

0.20 

3.0 

3 

0.66 

0.50 

0.30 

12.0 


Joseph, J. H: On the calculation of solar radiation fluxes in the troposphere. 
Solar Energy, Vol. 13, pp. 251-261. 


- 109 - 



XI— 5 


T 



The equivalent water vapor amount of the cloud layer (u* ) can be roughly 
determined by solving the following equation 

Aj -S 0 jj -A(u ct )J = (1 - Ri) S 0 [A(u ct + 1 .66 u c *) - A(u Ct )j 

Assumed values of the surface albedo (a s ) are 

bare land ... 0.14 sea . . . 0.07 snow ... 0.7 ice ... 0.4. 
The albedo of multi-layer clouds is determined as follows: 
i) 2 layers of cloud: 

Let the albedos of the two clouds be R x and R 3 . We need to consider 
the multiple reflection between the two cloud layers. The total transmisivity, T 8/ 
is the sum of T IS , T ia , T lS ^, ; where 

T,.' = (I -R,XI -R.) , 

T ia " - 0 - R.JO - R a ) R, R a , 
ij" = (1 -R l XI-R a )R*R. 3 , 
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I II ill + 

T ia " * ia + T ia + 4a * * 1 * 

= (1 -y(i -R a )/(l R a > / 

R 1S = 1 - T l2 - (R x + R 8 - 2R, R a )/(l - R a Ra ) • 

K, ' 



. .] 


'12 


ii) 3 layers of cloud: 

Let the albedos of the three clouds be R x , R 3 and R 3 . The total albedo 
of the 3 layers of cloud can be estimated as if they were 2 cloud layers for which the 
albedos are R l3 and R 3 . Then 

R x 33 = ^ " ^133 

(l -R 13 )0 -R 3 ) 


= 1 - 


1 - R l3 R 3 


(1 -Ri)(1 -Ra)(1 ~R 3 ) 

“ 1 - (R 1 R 3 +R 3 R3 + R 3 R 1 ) + 2R 1 R a Ra 
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Clear sky 

For the absorbed part of the solar radiation, the downward fluxes, 
S M , at the level i are 

S a o = s o Cl -A((u t - u Q ) sec £ ) ] 

S a8 = s o [1 " A((u t - u a ) sec £ ) ] 

S a4 =S 0 [1 “ A((u T - u 4 ) sec^) ] 

S 8 s =s 0 [1 - A(u t secO] • 


Then, the absorption of solar radiation in the layer between i-1 and i +1 , ASj , 
becomes: 

AS X — S a0 - S a3 , 

AS 3 ~ S a2 - S a4 , 

AS 5 = S a4 - S a s 

The absorption, by the earth's surface, of solar radiation reaching the surface is 

absorbed part: S a6 

s 

scattered part: S 0 (1 - a 0 )/(l - a 0 a s ) , 

where the denominator is the correction factor due to the multiple-scattering 


0 


>ao 




Saa 






a 


6 7777777 


i ! 

/777777777777 
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between the atmosphere and the earth's surface. The total solar radiation 
absorbed by the earth ' s surface is 

s 6 = (i -a s) [s. e + sl jf^aJLy]. 

Cloudy sky 

The general form for the optical length of the layer between the top of 
the atmosphere and level i, with the cloud layers, D Tl , is 

D t0 = (u T - u 0 ) sec£ , 

D Ta = d to + 0 ” CL X )(u 0 - u 3 ) sec £ +1 .66*CL 1 *u c * , 

D t 4 = d t 2 + 0 “ CL X )(1 - CL a )*(u s - u 4 ) sec £ 

+ 1.66 [CL,-u c * +CL t -(l -CL 3 )'(u 8 -uJ] , 

D t s - D T4 +(1 - CL 3 )-u 4 f(l - CL, xi -CL,) sec? 

+ [(1 -CLJ-Cl^+O -CL.KLJ- 1.66 } 

+ CLg (i% + u C 3 )• 1 .66 . 
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The downward solar radiation in the absorbed part, at level i 
(i = 0,2 ,4, 6), can be expressed as 

S„o = $o Cl-A(D t0 )], 

s., = s‘(l - CR, )• [ 1 -A(D t ,)] , 

S. 4 =S* (1 -CRJO - CR S )* [ 1 -A(D T4 )] , 

s. e = So (i -ago - cr 3 )(i -cr 3 )[1 -a(d, 8 )], 

where CR, = CL, -R, , and i = 1 ,2 and 3. Then, the absorptions in the 
atmospheric layers with clouds are 

AS, = (1 - CR, ) S. 0 - S. 3 , 

AS 3 = (1 - CR a ) S. a - S. 4 , 

ASg = 0 -CR 3 )S, 4 -S. s . 


In the absorbed part, some of the solar radiation, which is reflected 
at the tops of cloud 2 and cloud 3 (R a S a3 and R 3 S b4 ), can reach the earth's 
surface, in addition to S aQ , after multiple reflections between the three layers of 
cloud. If these amounts are defined by S a6 and S a6 , which refer to R„ S a3 and 
R 3 S a4 / respectively, they can be formulated as follows: 


S a 6 


S 


II 

a® 


b . M 1 -*•»> 

SaaR * 1 -R^Raa 7 
R i* (1 " Rs) 

Sa4Rs 1 -R 13 R 3 * 
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The above expressions are applied when all three layers have complete cloud 
coverage . To obtain more general formulae, the cloudiness must be taken into 
consideration. That can easily be done by using CR t and CR t 3 in place of R t and 
R , respectively. Then, the solar radiation in the absorbed part which reaches 
the earth's surface becomes 


s„, -[s*+s 


a 8 


CR a CR Z (1 - CR 33 ) CR s CR ia (1 - CR 3 ) _, 

+ S 


1 - CR,* CR a3 


a4 1 -CR 1S -CR 3 J 1 -a s -CR l83 


CR U = (CR, + CR 3 -2CR 1 *CR J )/(1 -CRj-CRj) 


In the scattered part, the total albedo of the clouds can be expressed 
by replacing R t by Rj . As a more general formula, we use CRj 1 instead of Rj . 

CR; = CLj-Rj' (! = 1,2,3), 

where Rj is the albedo of cloud i for the scattered part. Then 

- DI _ . (1 -CR,')(1 -CR' S )(1 - CRU 

CR i* 3 “ 1 - (CR ' -CR^ +CR^-CRJ + CR^-CRJ) + 2CRj-CR^CR^ 


M 


where 


2 - M - (CRJ + CRg + CR 3 ) + CRJ -CR^ -CR^ , 


M = (1 -CR^Kl -CRgKI -CR 3 ) . 


We assume that the total albedo of the atmosphere with clouds for the 
scattered part is 

a c = 1 - (1 - CR* a3 )(l - do) . 
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The scattered part which reaches the earth's surface is 

1 - a. 


*C S 

Sss = 1 Sr 


1 - a c - a s 

Finally/ the total solar radiation at the earth's surface becomes 
S as + S$s‘ Therefore, the absorption at the earth's surface is 


Se = (1 - a s )(S a3 + S ss ) 

1 


=(, ‘ as) {r^r 

1 - a 


S ’^^183: 


CR x CR a (1 - CR 33 ) CR 3 -CR 13 (1 -CR 3 ) 
I S a6 +S sa — r- /*n " /-»n + S. 


13 1 - CR t • CR 


a4 1 - CR 3 .CR 


1 - a c a s 


<}. 
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